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13.  ABSTRACT  (MiMimum  200  **otdi) 


Conventional  frequency  wave  number  -  processing  of  array  data  was 
used  to  estimate  the  azimuths  and  phase  velocities  of  fundamental  mode 
Rayleigh  waves  in  the  geology  and  ambient  noise  setting  at  North  Haverhill 
New  Hampshire.  These  estimates  arc  compared  to  true  azimuths  between 
known  shot  to  array  site  locations  to  determine  the  horizontal  refraction  of  Rg 
waves  due  to  the  regional  geology  and  the  variation  of  azimuth  and  phase 
velocity  estimates  due  to  additive  noise. 

D/v  „Thfc  ana,ysis  was  based  on  calculations  of  the  spectral  density  function 
t  Kx.Ky.f),  a  measure  of  power  as  a  function  of  frequency  and  wave  number 
tor  propagating  waves.  This  study  measured  0.3  to  1.0  second  period  Rg  waves 
with  phase  velocities  ranging  from  2.5  to  2.9  km/sec  having  refraction 
raypaths  with  a  lateral  velocity  gradient,  y  .  ^  /  • 
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1  INTRODUCTION 


This  thesis  concentrates  on  the  processing  of  seismic  data  acquired  as  pari  of 

the  New  York-New  England  Experiment  (NY-NEX)  conducted  from  September  17.  1988 
to  September  29,  1988.  NY-NEX  was  a  joint  venture  of  the  United  States  Air  Force,  the 
United  States  Geological  Survey  (USGS)  and  the  Geological  Survey  of  Canada  (GSC'i  to 
measure  the  seismic  velocity  structure  of  the  earth's  crust  and  to  test  seismic  methods 

that  could  be  used  to  verify  nuclear  test  ban  treaties  with  the  Soviet  Union.  This 

thesis  is  appropriate  for  test  ban  treaty  research  because  the  geology  of  the 
experiment  area  ia  analogous  to  that  found  in  the  Soviet  Union. 

NY-NEX  was  a  long  range  seismic  refraction  profile  experiment  that  consisted 

of  23  chemical  explosions  (shots)  located  at  sites  in  southern  Ontario,  northern  New 
York,  Vermont,  New  Hampshire,  and  central  Maine.  Each  of  the  shots  was  buried  and 
consisted  of  cither  2000  or  3000  pounds  of  explosives  with  the  exception  of  the 
easternmost  shot  which  w-as  6000  pounds. 

The  seismic  analysis  described  in  this  thesis  is  focused  on  estimating  azimuths 
and  phase  velocities  of  crustal,  fundamental  mode  Rayleigh  waves  (Rg  waves) 

generated  by  the  NY-NEX  shots  in  a  frequency  range  of  0.6  to  .3.3  Hz.  Rg  wave 
azimuths  and  phase  vclocitcs  were  measured  by  an  array  of  seismic  sensors  and 
compared  to  the  azimuths  and  group  velocities  calculated  from  known  shot  and  site 
coordinates  (listed  in  Table  1-1).  Rg  waves  arc  short-period,  fundamental  mode, 

dispersive  Rayleigh  surface  waves  that  propagate  along  a  free  surface  with  an 

elliptical  particle  motion.  The  paritclc  displacement  u(x,y,z,t)  of  a  fundamental  mode 

plane  surface  wave  propagating  in  the  x-y  plane  is  mathematically  represented  in 
cartesian  coordinates  flj  by 

u(x,y ,/.,()  =  A(z)  exp[j(Kx  x  +  Kv  y  -  wt)]  ( 1  ) 

where  Kx  and  Ky  arc  the  wavenumbers  in  the  x  and  y  directions,  w  is  l he  angular 
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frequency  (w  =  2nf,  where  f  is  the  frequency),  t  is  time,  and  A(z)  is  the  amplitude  as  a 
function  of  depth  and  frequency.  Rg  wave  velocities  vary  according  to  the  elastic  constants 
and  densities  in  the  near  surface  layers.  For  shallow  shots  in  New  England,  Kafka  [21  has 
shown  that  Rg  waves  travel  within  a  group  velocity  range  of  1.9  to  3.4  km/sec  in  the  0.6 
to  3.3  Hz  band.  The  group  velocity  U(f)  is  related  to  the  wavenumber  K  and  the  phase 
velocity  c(f)  by 

U(f)  =  Ry(tl-tO)  =  Aw/fcK  =  dw/dK  (2) 

and 

c(f)  =  ^R/ftm-tn)  =  w/K  (3) 

where  R  is  the  distance  from  shot  (origin  time  tO)  to  site  (arrival  time  tl),  and  ar  is  the 
distance  from  sensor  n  (arrival  time  tn)  to  sensor  m  (arrival  time  tm).  The 
frequency-wavenumber  (F-K)  power  spectral  density  function  P(Kx,Ky,f),  computed  from 
group-velocity  windows  of  the  Weston  Observatory  (WO)  array  data,  was  used  to  determine 
the  directional  and  phase  velocity  characteristics  of  Rg  waves  generated  from  the  NY-NEX 
shots.  P(Kx,Ky,f)  measures  power  as  a  function  of  frequency  and  wavenumber  for 
propagating  waves  and  can  be  used  to  estimate  azimuth  and  phase  velocity.  Deviations  of 
the  wave  direction  from  the  great-circle  source  azimuth  is  explained  by  (a)  lateral 
heterogeneity  in  the  shallow  crust  along  the  path  between  shot  and  receiver  and  (b)  the 
effects  of  additive  noise  on  those  estimates. 
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The  Weston  Observatory  (WO)  NY-NEX  array  was  located  i: 
North  Haverhill,  NH  with  its  center  at  latitude  44°  04' 
45.2"  (44.07922°),  longitude  72°  00'  30.7"  (72.00853°). 
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2  REGIONAL  GEOLOGY 


2.0  INTRODUCTION 

The  WO  NY-NEX  array  was  located  at  the  Dean  Memorial  Airport  in  North 
Haverhill,  New  Hampshire.  Since  the  area  is  flat,  the  array  analysis  can  be  confined  to 
the  x-y  coordinates.  The  earth  below  the  array  is  made  up  of  a  thick  overburden  of  glacial 
drift  and  alluvium.  Van  Diver  [31  describes  the  bedrock  geology  in  the  area  as 
characteristic  of  metamorphosed  sedimentary  and  volcanic  rocks  resulting  from  igneous 
activity.  Figures  1A,  IB,  and  2  show  the  large  scale  folds,  faults,  and  contacts  in  the 
region  as  having  a  north-northeast  trend. 


2  1  FORMATIONS 

The  French  Pond  granitic  (fpgc)  intrusion  is  north  of  the  WO  NY-NEX  array  (see 
Figure  1A).  Several  other  formations  that  lie  between  the  NY-NEX  shots  and  the  array- 
site  are  the  Albee  (Oal),  Bethlehem  gneiss  (bg),  Gile  Mountain  (Og),  Meetinghouse  Slate 
(Om),  Waits  River  (Ow),  and  the  Ammonoosuc  Volcanics  (Oam).  Table  2-1  describes  the 
formations  and  their  thicknesses  (where  known)  of  rocks  between  the  NY-NEX  shots  aid 
the  WO  receivers.  The  cross-section  in  Figure  IB  also  shows  these  bed  thicknesses  aid 
relative  dips  in  the  region. 


Table  2-1  NY-NEX  VT  -  NH  ROCK  FORMATIONS 


F CRM  DESCRIPTION  THICK 

_  (km) 

eg  Granodiorite  or  quartz  monzonite 

fpgc  Porphyritic  to  coarse  grain  pink  biotite  granite 
Cal  Feldspathic  quartzite  in  chlorite  zone  1.5 

0  a  m  S  oda-rhyolite  and  chlorite  schist  2.5 

Zg  Gradational  mica  schist  to  quartz-mica  schist  2.0 


Thinly  bedded  mica  and  quartz-mica  slate  <&  schist  0.8 
Interbedded  calcareous  and  non-calcareous  schists  3.0 


2.2  GEOLOGICAL  HISTORY 

The  geologic  history  of  surface  rocks  in  the  study  area  began  in  the  Middle 
Ordovician  with  deposition  of  the  Waits  River,  Gile  Mountain.  Meetinghouse,  and  Alive 
formations  [4|.  Minor  deformation  of  these  rocks  resulted  from  the  Taconic  orogeny. 
This  deformation  is  noticeable  in  areas  such  as  Littleton,  NH,  where  Silurian  rocks  rest 
uncomformabiy  on  top  pre-Silurian  rocks.  Major  tectonic  deformation  during  the  Middle 
Devonian  occurred  in  the  Acadian  Orogeny,  resulting  in  the  folding,  faulting, 
metamorphism,  and  numerous  igneoLi  intrusions  such  as  the  French  Pond  granite  (fpgc). 
The  middle  Ordovician  Meetinghouse  slate  and  Gile  Mountain  formations  of  Vermont  are 
separated  from  the  Silurian  to  middle  Devonian  Albee  formation  of  New  Hampshire  by  the 
north-northeast  trending  Monroe  fault  (41.  East  of  the  Monroe  fault  lies  the  Ammonoosuc 
thrust  fault  that  dips  36°W.  The  Ammonoosuc  thrust  fault  establishes  the  contact  between 
the  New  Hampshire  Albee  formation  and  the  Ammonoosuc  Volcanics. 

Recent  geologic  changes  are  attributed  to  glaciation  and  stream  erosion.  Wisconsin 
glaciation  (18,000  to  20,000  years  ago)  has  the  only  lasting  glacial  effect  on  the 
geomorphology.  This  glacial  erosion  is  seen  in  the  rounded  valleys  and  carved  flanks  of 
the  Presidential  Range.  Stream  erosion  is  mostly  attributed  to  the  Connecticut  River  which 
winds  through  the  Ordovician  Albee  formation  (Oal). 

2.3  GEOPHYSICAL  ROCK  PROPERTIES 

The  regional  Rayleigh  wave  velocity  estimates  (Vr)  for  rocks  in  the  formations 
located  between  NY-NEX  shot  and  WO  receiver  points  are  listed  in  Table  2-2.  Rayleigh 
wave  velocities  (Vr),  listed  in  Table  2-2,  were  calculated  from  the  equations 

Vr  =  0  9194  Vs  (4) 

and 
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Vs  =  Vp/V3 


(5) 


where  Vp  and  Vs  are  the  compressional  and  shear  wave  velocities.  Equations  (4)  and 
(5)  are  based  on  a  Poisson’s  ratio  of  0.25  and  are  further  discussed  by  Sheriff  et  al.  [5], 
These  estimates  of  Vr  provide  a  velocity  comparison  between  rock  formations.  Their 
calculations  are  based  a  model  where  the  travel  path  is  along  the  surface  (zero  depth)  of 
a  homogeneous  medium  (which  does  not  exhibit  dispersion  [5]).  Values  for  Vp  and  Vs 
were  found  in  the  Handbook  of  Physical  Constants  [6],  The  handbook  estimates  a  quality 
factor  (Q)  for  rocks  in  the  region  ranging  between  125  and  314,  which  is  characteristic 
of  the  higher  attenuation  levels  found  in  igneous  and  metamorphic  rocks. 

To  find  the  relationship  between  the  rock  formations  and  the  rock  velocities 
surrounding  the  WO  array,  seismic  velocities  for  the  bg,  fpgc,  Oal,  Oam  and  Ow 
formations  were  estimated  using  Tables  2-1  and  2-2.  The  first  significant  formation  is 
so^.neast  of  the  array  and  is  bg;  it  is  composed  of  granodiroite  and  is  assumed  to  have 
a  velocity  Vr  of  2.3  km/sec.  The  fpgc  intrusion  to  the  north  of  the  array  mostly  consists 
of  granite  and  therefore  is  assumed  to  have  a  velocity  Vr  of  2.7  km/sec.  The  Og 
formation  is  dominated  by  schist  and  quartz-mica  schist  and  is  thought  to  have  a  velocity 
Vr  of  3.0  km/sec.  Oal  is  assumed  to  have  a  velocity  Vr  of  about  3.0  km/sec  since  it  is 
composed  of  the  quartzite  and  chlorite  schist  rocks  listed  in  Table  2-2.  Oam  is  in  the 
chlorite  zone  and  expected  to  be  a  low  velocity  layer  with  a  Vr  of  about  2.5  km/sec. 

Rg  wave  group  velocities  in  the  NY-NEX  study  area  are  therefore  expected  to  be 
in  the  range  of  2.5  to  3.0  km/sec,  characteristic  of  the  high  velocity  metamorphic  and 
igneous  formations  in  the  area  For  example,  Figure  2  shows  the  geologic  structure 
between  shot  7  (West  Groton,  Vt)  and  the  WO  array  consisting  of  a  cross-section  of  the 
Ow,  Og,  Om,  and  Oal  formations.  The  Rg  rock  velocities  for  these  formations  are 
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expected  to  be  in  the  range  of  2.5  to  3.0  km/sec. 


Table  2-2  ROCK  VELOCITIES  AT  10  BARS  PRESSURE 


Rock 

Vp 

Vs 

Vr 

Measure  Area 

(km/sec) 

(km/sec) 

(km/sec) 

Alluvium 

0. 5-2.0 

0.3-1. 2 

0.3-1 .1 

Glacial  till 

0. 4-1.0 

0.2-0. 6 

0 . 1-0 . 6 

Granite 

5.1 

2.9 

2.7 

Barre, Vt 

Granodiorite, gneiss  4.4 

2.5 

2.3 

Bethlehem, NH 

Chlorite  Schist 

4  .  8 

2.8 

2.5 

Chester 

Quartzite 

_ 

4  .  0 

3.7 

Quarry  Vt 

Slate 

5 . 5 

3.2 

2.9 

Medford, Ma 

Vr  is  the  Rayleigh  wave  velocity  along  the  surface 
(zero  depth)  of  an  isotropic,  homogeneous  medium. 


3  DATA  ACQUISITION 


3.0  INTRODUCTION 

Shot  times  for  the  NY-NEX  experiment  were  restricted  to  the  early  morning  hours 
to  take  advantage  of  the  low  background  seismic  noise  conditions  resulting  from  lulls  in 
cultural  activity  and  benign  weather  patterns.  The  seismic  profile  encompassed  three  nights 
of  shot  sequences.  These  shot  sequences  were  separated  on  each  night  by  a  1.0  to  1.5 
hour  interval.  Each  shot  sequence  consisted  of  four  or  five  shots  at  different  locations 
exploded  at  two  or  three  minute  intervals.  Each  night  of  shooting  was  followed  by  a  three 
to  four  day  rest  period.  Table  3-1  lists  the  NY-NEX  shot  sizes,  shot  times,  and  shot 
locations.  A  geographic  map  of  the  shot  locations  is  shown  in  Figure  3. 

All  of  the  shot  sequences  were  digitally  recorded  on  magnetic  tape  and  transferred 
to  WO  where  they  were  divided  into  two  minute  windows  for  data  analysis.  Each  window 
begins  at  the  shot  time.  A  pre-shot  window  containing  ambient  seismic  noise  was  also 
recorded  prior  to  each  shot  sequence. 

3  1  SITE  INFORMATION 

The  WO  NY-NEX  array  was  located  approximately  2  km  southeast  of  North 
Haverhill,  New  Hampshire  at  the  Dean  Memorial  Airport  (see  Figure  4).  The  geographic 
coordinates  of  the  center  of  the  array  was  located  at  latitude  44.07923°  (44°  04’  45.2"N) 
and  longitude  72.00852°  (72°  00’  30.7"W). 

The  variably  spaced  NY-NEX  array  was  designed  to  suppress  ambient  seismic 
noise  and  to  capture  a  wide  velocity  range  of  NY-NEX  seismic  signals.  The  shape  of 
the  array  represents  a  compromise  between  maintaining  linear  legs  in  the  array  and  using 
available  airport  grounds  for  sensor  planting.  Two  array  configurations  were  used  to  record 
NY-NEX  signals.  Sensor  coordinates  for  these  configurations  are  listed  in  Table  3-2  and 
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shown  in  Figures  5A  and  5B.  The  only  difference  between  the  two  arrays  is  that  in  the 
"short”  array  (Figure  5B)  sensor  #1  was  moved  106.6  meters  (m)  closer  to  sensor  #2.  This 
change  was  made  to  reduce  wind  noise  generated  by  the  line  of  trees  next  to  the  sensor 
ai  the  north  end  of  the  airfield.  Most  of  the  seismic  sensors  were  buried  on  level  ground 
at  173.7  m  above  sea  level  with  the  exception  of  channels  14  and  15  which  were  buried 
at  elevations  of  170.7  m  and  164.6  m  above  sea  level  as  listed  in  Table  3-2  and  shown 
in  Figure  4.  The  center  of  the  array  was  located  at  the  3-component  station  of  /'hannels 
7,  8  and  9.  Channels  1-7  and  10-15  define  the  two  main  legs  of  the  array.  Channels  1-7 
comprised  the  N  8°  1’  W  leg  of  the  array.  Channels  10-15  comprised  the  N  69°  26.5’  W 
leg  of  the  array.  Channel  16  bisected  the  angle  between  the  two  legs  and  was  69.5  m 
from  the  center  of  the  array  (see  Figures  5 A  and  5B).  The  placement  of  sensor  16  was 
based  on  the  analysis  of  a  series  of  spatial  window  function  models.  The  spatial  window 
function  (see  Section  3.2.3)  for  the  final  array  pattern  is  shown  in  Figure  6. 
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Table  3-1  SHOT  INFORMATION 


JULIAN 

DAY 

TIME 

(HHMM) 

SHOT 

SIZE 

(lbs) 

SHOT  LOCATION 

261 

0400 

2 

2000 

Farmington  Falls,  ME 

261 

0402 

5 

2000 

Jefferson  Notch,  NH 

261 

0404 

7 

3000 

West  Groton,  NH 

261 

0406 

22 

2000 

Warner,  NY 

261 

0408 

14 

3000 

Long  Lake,  NY 

261 

0600 

6 

2000 

Littleton,  NH 

261 

0602 

4 

2000 

Gilead,  ME 

261 

0604 

1 

6000 

Waterville,  ME 

261 

0800 

3 

2000 

Mount  Zircon,  ME 

261 

0802 

23 

2000 

Sanford,  ME 

261 

0804 

10 

3000 

West  Addison,  VT 

268 

0400 

8 

2000 

Barre  West,  VT 

268 

0402 

9 

2000 

Lincorn  Gap,  VT 

268 

0404 

12 

2000 

Newcomb,  NY 

268 

0406 

22 

2000 

(repeat) 

268 

0408 

20 

3000 

Marmora,  Ontario 

268 

0600 

7 

3000 

(repeat) 

268 

0602 

17 

3000 

Alexandria  Bay,  NY 

268 

0604 

13 

2000 

Newcomb,  NY 

268 

0606 

10 

2000 

(repeat) 

268 

0800 

14 

3000 

(repeat) 

268 

0802 

11 

2000 

North  Hudson,  NY 

268 

0804 

21 

2000 

Londonderry,  VT 

268 

0806 

4 

3000 

(repeat) 

274 

0400 

20 

2000 

(repeat) 

274 

0402 

18 

2000 

Sharpton,  Ontario 

274 

0404 

17 

500 

(repeat) 

274 

0406 

14 

3000 

(repeat) 

274 

0600 

19 

2000 

Lime  Lake,  Ontario 

274 

0602 

16 

2000 

Fort  Drum,  NY 

274 

0604 

15 

2000 

Star  Lake,  NY 

274 

0606 

10 

3000 

(repeat) 

All  shots  occurred  within  0.015  seconds  of  the  listed 
Greenwich  Mean  Time  (GMT) . 
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Table  3-2  SENSOR  COORDINATES 


CHANNEL 

SERIAL 

NUMBER 

X-LOC 

(m) 

Y-LOC 

(m) 

Z-LOC 

(m) 

RANGE 

(m) 

1 

363Z 

-62.5 

443.6 

173.7 

448 . 0 

2 

347Z 

-33.7 

239.0 

173.7 

241.4 

3 

349Z 

-18 . 1 

128.2 

173.7 

129.5 

4 

366Z 

-9.7 

68.8 

173.7 

69 . 5 

5 

279Z 

-5.2 

36.9 

173.7 

37.3 

6 

322Z 

-2.8 

19.8 

173 . 7 

20.0 

7 

348Z 

0.0 

0.0 

173.7 

0.0 

8 

257N 

0.0 

0.0 

173.7 

0.0 

9 

258E 

0.0 

0.0 

173.7 

0.0 

10 

278Z 

-18.7 

7.0 

173.7 

20.0 

11 

367Z 

-34.9 

13.1 

173.7 

37.3 

12 

356Z 

-65.1 

24.4 

173.7 

69.5 

13 

32 1Z 

-121.3 

45.5 

173.7 

129.5 

14 

319Z 

-226.0 

84.8 

170.7 

241.4 

15 

346Z 

-391.7 

119.9 

164 . 6 

341.4 

16 

362Z 

-43.9 

53.8 

173.7 

69.5 

Position  change  of  channel  1  on  day  267  at  1900  hours 

CHANNEL  SN  X-LOC  Y-LOC  Z-LOC  RANGE 

1  363Z  -47.6  338.1  173.7  341.4 

The  NY-NEX  array  center  (channel  7,8,9)  was  located 
at  latitude  44°  4'  45.2"  (44.07922°)  and  longitude 
12°  00'  30.7"  (72.00853°).  X-LOC,  Y-LOC,  and  Z-LOC 
represent  the  distance  north,  east,  and  depth  of  the 
sensor  positions  from  the  array  center. 
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3.2  RECORDING  HARDWARE 


3.2.1  Geophysical  Data  Acquisition  System 

The  Geophysical  Data  Acquisition  System  (GDAS)  developed  at  WO  provided  the 
hardware  [7]  and  software  [8)  to  record,  store,  and  process  16  channels  of  data  from  the 
NY-NEX  experiment.  GDAS  is  a  stand-alone  system  based  on  Digital  Equipment 
Corporation’s  Q-busT  microcomputer.  The  18-bit  Q-bus  supports  the  LS1-11/231 
microprocessor  used  for  data  processing  and  interrupt-driven  data  acquisition.  Digitized 
data  were  stored  on  floppy  disks  and  9-track  magnetic  tapes  for  further  (in-house)  data 
analysis.  A  15 -bit  analog-to-digital  converter  (A/D)  digitizes  analog  filtered  sensor  data. 
A  diagram  of  the  GDAS  computer  is  presented  in  Figure  7.  For  the  NY-NEX  experiment, 
a  niter  card  for  each  channel  consisted  of  a  6-pole  low-pass  uutterwunh  filter  at  34.3  Hz 
and  a  2-pole  low-pass  butterworth  filter  at  100.0  Hz.  Filter  cards  contained  a  hardware 
modified  gain  that  was  set  at  a  high  or  low  level  according  to  the  shot  to  receiver  distance 
(see  Table  3-3).  The  roll-off  of  the  butterworth  filter  for  the  low  and  high  gain  filter  cards 
is  presented  in  the  system  transfer  function  (channel  response)  plots  shown  in  Figure  8. 
The  GDAS  field  unit  was  hardware  switch  selected  to  record  data  at  a  100  Hz  sample  rate. 

The  WO  NY-NEX  array  deployed  2  horizontal  and  14  vertical  Electro-tech  EV-17 
short-period  (1  Hz)  seismometers  featuring  an  inclined  spring  suspension  with  a  moving 
magnet  mass  of  2760  g.  The  seismometers  have  a  nominal  sensitivity  of  8.8  volts/mm/sec. 
Ground  motions  were  converted  by  the  seismometers  into  voltages  which  were  transmitted 
to  the  preamplifiers.  These  preamplifiers  conditioned  the  signals  to  desired  levels  and 
transmitted  them  to  the  analog  filter  cards. 


T  Q-bus  and  LSI- 11/23  are  trademarks  of  Digital  Equipment  Corporation. 
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Table  3 

z2 

SHOT  RECORDS 

SHOT 

SIZE 

(lb) 

TIME 

(GMT 

} 

FILE 

NAME 

GAIN 

ARRAY 

SIZE 

SHOT 

DATE 

3 

2000 

261 

08 

00 

W3D61B 

HIGH 

LONG 

17-SEP-83 

4  A 

2000 

261 

06 

02 

W2D61C 

LOW 

LONG 

17-SEP-88 

4B 

3000 

268 

08 

06 

W3D68E 

HIGH 

SHORT 

24-SEP-88 

5 

2000 

261 

04 

02 

W1D61C 

LOW 

LONG 

17-SEP-38 

6 

2000 

261 

06 

00 

W2D61B 

LOW 

LONG 

17-SEP-88 

7A 

3000 

261 

04 

04 

W1D61D 

LOW 

LONG 

17-SEP-88 

7B 

3000 

268 

06 

00 

W2D68B 

LOW 

SHORT 

24-SEP-88 

8 

2000 

268 

04 

00 

W 1 D  6  8  B 

HIGH 

SHORT 

24-SEP-88 

9 

2000 

268 

04 

02 

W1D68C 

HIGH 

SHORT 

24-SEP-38 

10A 

3000 

261 

08 

04 

W3D61D 

HIGH 

LONG 

17-SEP-88 

10B 

2000 

268 

06 

06 

W2D68E 

LOW 

SHORT 

24-SEP-88 

IOC 

3000 

274 

06 

06 

W2D74E 

LOW 

SHORT 

29-SEP-88 

11 

2000 

268 

08 

02 

W3D68C 

HIGH 

SHORT 

24-SEP-85 

12 

2000 

268 

04 

04 

W1D68D 

HIGH 

SHORT 

24-SEP-88 

13 

2000 

268 

06 

04 

W2D68D 

LOW 

SHORT 

24-SEP-88 

14A 

3000 

261 

04 

08 

W1D61F 

LOW 

LONG 

17-SEP-88 

14B 

3000 

268 

08 

00 

W3D68B 

HIGH 

SHORT 

24-SEP-38 

14C 

3000 

274 

04 

06 

W1D74E 

HIGH 

SHORT 

29-SEP-88 

21 

2000 

268 

08 

04 

W3D68D 

HIGH 

SHORT 

24-SEP-88 

22A 

2000 

261 

04 

06 

W1D61E 

LOW 

LONG 

17-SEP-88 

22B 

2000 

268 

04 

06 

W1D68E 

HIGH 

SHORT 

24-SEP-88 

23 

3000 

261 

08 

02 

W3D61C 

HIGH 

LONG 

17-SEP-88 

A  table  of  data  . 
array  site  range 

analyzed  within  a 
of  200  km. 

NY-NEX 

shot  to  WO 
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The  filtered  outputs  were  transmitted  to  the  A/D  where  they  were  multiplexed,  digitized, 
and  stored.  This  process  is  illustrated  in  the  following  diagram  for  each  of  the  digitized 
signals  by 

g(t)  +  n(t)  ->  SENSOR  ->  PREAMPLIFIER  ->  FILTER  ->  A/D  ->  s(nt) 
where  g(t)  +  n(t)  represents  the  analog  recorded  signal  plus  noise  at  the  sensor  and  s(nt) 
is  the  digitized  recorded  signal  at  the  sampling  interval  nt.  The  sampling  sequence  of  the 
signal  is  represented  mathematically  by 

s(nt)  =  p(nt)  (g(t)  *  h(t)l  (6) 

where  *  is  the  symbol  for  convolution,  p(nt)  is  defined  as  the  digitizer  or  the  Dirac  comb 
function,  and  h(t)  is  the  time  domain  representation  of  the  channel  response  (further 
discussed  in  Section  3.3.2).  Typical  examples  of  both  raw  and  filtered  ground  motions 
g(nt)  +  n(nt)  recorded  by  the  NY-NEX  array  are  shown  in  Figure  9. 

3.2.2  Channel  Calibration 

Channel  calibration  procedures  are  used  to  describe  the  system  mathematically  in 
terms  of  a  system  transfer  function.  In  this  thesis,  the  system  transfer  function  is  identified 
for  each  channel  as  the  combined  measurement  of  seismometer,  preamplifier,  and  analog 
filter  card  electronics.  Since  the  system  transfer  function  includes  the  sensor  electronics, 
it  is  herein  referred  to  as  the  channel  response. 

The  first  step  in  system  calibration  is  the  application  of  a  step  in  electric  current 
(generated  by  the  digital  to  analog  convener)  to  the  seismometer  calibration  coil.  This 
causes  an  equivalent  weight  lift.  The  output  (calibration)  wavelet  due  to  this  "weight  lift" 
is  recorded  by  GDAS.  This  procedure  is  represented  mathematically  for  the  seismometer 
in  terms  of  equivalent  force  F  by 

F  =  Ma  =  mg  =  Gi  (7) 
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where  M  is  the  mass  of  the  seismometer  weight,  a  is  the  acceleration  of  the  seismometer 
weight,  m  is  the  calibration  mass,  g  is  gravity  (9.81  m/sec2),  G  is  the  motor  constant,  and 
i  is  the  calibration  current. 

The  second  step  in  system  calibration  is  to  match  the  calibration  wavelet  with  a 
computer  generated  second  order  seismometer  model  (of  known  electronics)  with  a  gain, 
phase,  and  damping  coefficient  [8].  A  match  is  selected  when  the  sum  of  squared  errors 
(SSE)  between  the  model  wavelet  and  the  calibration  wavelet  is  less  than  0.001,  or  when 
the  ambient  noise  in  the  calibration  wavelet  is  less  than  the  residual  (SSE  between  the 
model  and  calibration  wavelet)  value. 

Calibration  measurements  were  performed  on  site  during  GDAS  operation  for  both 
the  low  and  high  gain  filter  card  electronics.  To  improve  channel  calibration  wavelet 
estimates,  several  calibration  tests  were  recorded  and  averaged  for  the  low  and  high  gain 
electronics.  The  low  and  high  gain  channel  responses  are  listed  in  Table  3-4  and  shown 
in  Figure  8  for  the  NY-NEX  site.  The  amplitude  response  is  constant  to  ground  panicle 
velocity  in  the  2.0  to  30.0  Hz  range.  At  frequencies  greater  than  30.0  Hz,  there  is  a  sharp 
fall-off  in  the  channel  response  due  to  the  analog  filter. 
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Table  3 

-4  LOW  GAIN 

CHANNEL  PARAMETERS 

CHANNEL 

POLE  FREQUENCY 

DAMPING 

MODEL 

(Hz) 

(%) 

Amplitude  Ratio 

1 

0 . 932129 

0.690000 

-68.0309 

2 

0 . 969772 

0 . 648438 

-68 .4867 

3 

1 .010064 

0 . 613281 

-57.2372 

4 

0 . 937192 

0 . 669531 

-66 . 0749 

5 

0 . 956557 

0.71875 

-76.4375 

6 

0 . 947685 

0.683594 

-71  .2204 

7 

0 . 943213 

0 . 739844 

-78 .0973 

6 

0 . 996625 

0.641406 

-67 .0136 

9 

0 . 977805 

0.641406 

-66 . 6782 

10 

0 . 97739 

0 . 627344 

-65 . 652 

11 

0 . 958049 

0.704688 

-65.464 

12 

0 . 939770 

0.725781 

-75 . 6339 

13 

0 . 939239 

0.71875 

-76.5301 

14 

0 . 926679 

0 .746875 

-80.7887 

15 

0 . 926619 

0 . 697656 

-71  .  6362 

16 

0 . 929452 

0.6625 

-68 .1599 

HIGH  GAIN 

CHANNEL  PARAMETERS 

CHANNEL 

POLE  FREQUENCY 

DAMPING 

MODEL 

(Hz) 

<%> 

Amplitude  Ratio 

1 

0 . 937668 

0.693262 

-141.4317 

2 

0.969542 

0.640527 

-140 . 8091 

3 

1.009143 

0 . 610000 

-119.7630 

4 

0 . 934152 

0 . 657227 

-136.7298 

5 

0 . 956689 

0.718750 

-159.2683 

6 

0 . 940968 

0 . 662500 

-146 . 9642 

7 

0.935248 

0.725000 

-162 . 6823 

8 

1 . 008804 

0 . 655469 

-141  .  8040 

9 

0 . 985984 

0.662500 

-141  .4921 

10 

0 . 974676 

0.606250 

-134.2870 

11 

0 . 953672 

0.698535 

-135 . 1754 

12 

0 . 938281 

0.720000 

-156 .7691 

13 

0 . 936970 

0 . 695000 

-156.4051 

14 

0 . 931250 

0.750000 

-168.0947 

15 

0 . 927178 

0 . 687500 

-148.5644 

16 

0.924888 

0.649316 

-140 . 0599 
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3.2.3  Array  Attributes  (Array  Calibration) 


The  spatial  window  function  I  B(K)  1 2,  also  referred  to  as  the  beamforming  array 
response,  was  calculated  as  the  magnitude  of  the  spatial  transform  of  the  array  pattern's 
response  to  a  unit  impulse  of  infinite  phase  velocity.  A  continuous  array  response  is 
calculated  by 

oo  oo 

B(Kx,Ky;0  =  11  w(x,y)  exp(-j(Kx  x  +  Ky  y)j  dx  dy  (8) 

— OO  -oo 

Equation  (8)  for  an  array  of  NS  sensors  with  wavenumbers  Kx  and  Ky  is  digitally 
represented  at  sensor  locations  m  and  n  by 

NS 

B(Kx,Kv;t)  =  1/NS:  I  Wmn  exp{-j[Kx(xn-xm)  +  Ky(yn-ym)])  (9) 
m,n=l 

The  weight  function  Wmn  in  Equation  (9)  allows  weighting  for  array  optimization 
[9,10.11,12!  Weighting  is  itself  a  unique  venture  of  finding  the  correct  sensor  weighting 
to  optimize  azimuth  and  phase  velocity  estimation  and  is  not  exploited  in  this  thesis.  For 
a  uniformly  weighted  function  of  Wmn=l,  the  spatial  window  function  lB(K)l:  is 
calculated  using 

NS 

B(Kx,Ky,t)  =  I  1/NS  I  exp[-j(Kx  xn  +  Ky  yn))  1 2  (10) 

n=l 

Figures  6A  and  6B  show  the  uniformly  weighted  spaiial  window  functions  for  the 
two  NY-NEX  array  patterns.  The  data  are  displayed  over  a  wavenumber  range  of  ±0.1 
rad/m  with  contour  steps  in  decibels  (db)  calculated  by  -101og[  lB(K)l2],  where  I  B(K)  1 2 
has  been  normalized  by  its  maximum.  Table  3-5  provides  insight  into  the  relation  between 
db  and  the  percentage  of  the  maximum  normalized  beamforming  array  response. 

The  beams  shown  in  Figures  6A  and  6B  will  optimally  detect  waves  from  the  NW 
and  SW  directions.  Rg  wave  azimuth  and  phase  velocity  are  shown  in  Section  5.1.2  to 
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be  calculated  from  the  Kx  and  Ky  values  where  P(Kx,Ky,0  is  a  maximum. 


Table  3-5 

DECIBEL 

CONTOUR  LEVELS 

Decibels 

%  of  Maximum 

0 .00 

100.00 

0.25 

94.41 

0.50 

89.13 

1.00 

79.43 

2.00 

63.10 

3 . 00 

50 . 12 

4.00 

39 . 81 

5.00 

31 . 62 

6 .00 

25.12 

7.00 

19.95 
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4  NOISE 


4.0  INTRODUCTION 

All  unwanted  recorded  information  is  defined  here  as  noise.  Noise  includes 
digitization  error,  aliasing,  GDAS  electronic  noise,  ambient  seismic  motion  (microseisms), 
scattered  refractions,  and  in  this  thesis,  all  that  is  not  the  Rg  wave. 

Ambient  ground  motion  may  be  defined  as  the  sum  of  independent  random  seismic 
events.  Random  events  may  be  described  as  noise  generated  by  vehicles,  animals,  and  all 
random  electronic  noise  not  described  by  the  channel  response.  Typical  noise  levels 
recorded  during  pre-shot  windows  are  shown  in  Figure  10.  Noise  can  be  described  both 
at  a  point  (seismometer)  and  that  which  Is  passed  by  the  array.  Random  noise  is  classified 
by  its  stationarity,  normality,  and  spectra.  Microseisniic  activity  is  found  to  exhibit  the 
attributes  of  a  stationary  Gaussian  random  process. 

Aki  and  Richards  [1]  state  that  the  scatter  in  phase  velocity  measurement  is  due 
to  lateral  heterogeneity  and  ambient  seismic  noise.  Aki  and  Richards  [1]  define  the  root- 
mean-square  (RMS)  phase  velocity  measurement  error  between  two  sensors  that  share  the 
same  signal  and  noise  as 

Acrms/c  =  l/(2n)lN(f)l  /  IS(f)l  (f>x)  (11) 

and  the  corresponding  azimuth  measurement  error  as 

«Drms  =  l/x/2  lN(f)l  /  IS(f)l  (12) 

where  I  N(f)  I  is  the  amplitude  spectra  of  the  noise,  IS(f)l  is  the  amplitude  spectra  of  the 
signal,  £  is  the  wavelength,  and  *x  is  the  distance  between  seismic  stations. 

4.1  NOISE  TEST 

The  Gaussian  probability  distribution  function  was  used  to  test  the  normality  of 
the  data  in  the  NY-NEX  pre-shot  ambient  noise  windows.  For  noise  to  be  classified  as 
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stationary,  it  must  be  independent  of  time.  The  Sax  test  [13]  for  a  stationary  normal 
process  was  used  to  test  the  behavior  of  the  ambient  noise  data. 

4,1.1  Test  for  Normality 

The  central  limit  theorem  states  that  the  sum  of  a  number  of  random  independent 
variables  has  a  distribution  which  approaches  a  Gaussian  function  as  the  number  of 
variables  becomes  large.  Mathematically,  a  random  variable  x  described  as  a  sum  of 
random  independent  events  ri 

x  =  rl  +  r2  +  r3  +  r4  +  ...  m  (13) 

will  approach  a  Gaussian  distribution  function  as  the  number  of  contributions  becomes 
large  [14].  Its  probability  density  function  p(x)  is 

p(x)  =  l/[cW(27t)]  exp[-(x-jx)2/(2crz)]  (14) 

where  |t  is  the  mean  of  x  and  o  is  the  standard  deviation.  Distributions  of  pre-shot 
ambient  seismic  motion  were  plotted  on  probability  paper  (as  shown  in  Figure  12)  to 
determine  if  the  noise  was  distributed  as  a  normal  Gaussian  process.  The  advantage  of 
the  probability  paper  is  that  it  transforms  the  distribution  of  the  Gaussian  variate  x’  into 
a  straight  line  [15],  and  the  test  becomes  one  of  accepting  a  linear  fit.  The  test  data  are 
based  on  calculating  the  variate  x’  by  inversely  solving  Equation  (15) 

x’ 

y  =  P(x’)  =  Prob(x  <  x’)  =  \  p(x)  dx  (15) 

-OO 

for  x’, 

x’  =  P'(y)  (16) 

P(x’)  is  the  probability  that  the  random  variable  x  is  less  than  or  equal  to  the  variate  x’. 
The  distributions  were  computed  by  dividing  observations  in  the  ambient  noise  window 
into  128  equal  histogram  bins  (class  intervals). 
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4.1.2  Test  for  Stationaritv 


For  stationary  Gaussian  noise,  spectral  coefficients  in  a  narrow  band  are  distributed 
as  a  Rayleigh  density  function  (exponential  probability  function) 

p(x)  =  x/( c2)  exp[-x2/(2c2))  for  x  >  0  (17) 

where  c  is  a  constant  such  as  the  standard  deviation  a.  Sax  [131  showed  that  the 
probability  of  a  Power  Spectral  Density  (PSD)  measurement  in  a  band  E  of  Gaussian  noise 
between  E  and  E+dE  is  described  by  the  exponential 

P(E)dE  =  exp(-E/E)  dE  (18) 

which  follows  the  Rayleigh  distribution  function 
E  E 

I  P(E)dE  =  \  exp[-x7(2c2))  (19) 

-oo  -oo 

A  linear  trend  in  the  graph  identifies  the  noise  as  stationary  and  Gaussian.  For 
the  test,  spectral  coefficients  were  generated  at  frequencies  in  the  0.6  to  3.3  Hz  band  (the 
band  where  the  Rg  wave  is  known  to  occur).  These  coefficients,  or  observations,  were 
plotted  as  linear  energy  (volts2/Hz)  versus  logarithmic  frequency  of  occurrence  (transformed 
scores)  as  shown  in  Figure  11.  The  spectral  coefficients  were  generated  by  calculating  34 
fast  Fourier  transforms  (FFT’s)  of  128  individual  samples  (plus  128  zero  pads  -  see 
Appendix  A)  in  the  same  ambient  seismic  noise  window  used  in  Table  4-1.  Spectral 
coefficients  were  gathered  at  frequencies  in  the  0.6  to  3.3  Hz  band  and  distributed  as  a 
cumulative  process. 
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4.1.3  Summary  of  the  Test  Results 


Table  4-1  summarizes  the  Gaussian  distribution  test  calculated  for  four  individual 
pre-shot  ambient  noise  windows.  Table  4-1  identifies  some  sensors  within  the  array  as 
non-Gaussian.  The  non-Gaussian  distributions  of  channels  14  and  15  are  probably  the 
result  of  these  channels  being  close  to  the  trees.  The  non-Gaussian  distributions  are 
probably  due  to  an  increase  in  the  wind  level,  as  seen  when  comparing  the  amplitude 
levels  of  channels  3  and  14  in  Figure  12. 


Table  4-1  GAUSSIAN  DISTRIBUTION  TEST 


CHANNELS 


FILE 

I 

2 

3 

4 

5 

6 

7 

10 

11 

12 

13  14 

15 

16 

W1D61A 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

W2D61A 

X 

X 

X 

X 

X 

0 

X 

X 

X 

0 

W2D68A 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

o 

W2D74A 

X 

X 

X 

X 

X 

X 

X 

o 

X 

X 

X 

X 

X 

X 

An  "x 

"  marks 

the 

channels 

that 

fit 

a 

Gaussian 

distribution . 

An 

"o 

"  marks  the 

:  channels 

that  are 

very 

close 

to 

a 

Gaussian 

distribution . 

A  " 

II 

(blank) 

marks 

the  channels  as  non-Gaussian.  File  W1D61A  was  based  on 
a  time  window  from  35.84  to  122.87  sec.  Files  W2D61A, 
W2D68A,  and  W2D74A  were  based  on  time  windows  from 
0.01  to  62.00  sec. 


Examples  of  the  stationarity  test  for  channel  3  are  shown  in  Figure  11.  The  plots 
show  that  the  channel  tests  as  non-stationary  at  1.76  and  2.15  Hz.  Note  that  these 
frequencies  also  exhibit  higher  amplitude  levels,  indicating  that  the  non-linearity  in  the 
plot  is  probably  due  to  wind  noise. 
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4.2  CHANNEL  NOISE 


4.2. 1  Hardware  Noise 

Channel  hardware  noise  provides  a  lower  bound  estimate  of  the  overall  system 
noise.  Channel  hardware  noise  consists  of  the  combined  measurement  of  electronic, 
seismometer,  and  digital  noise.  Hardware  noise  is  isolated  by  calculating  the  temporal 
incoherence  between  collocated  sensors.  The  noise  Gnn(f)  is  calculated  in  terms  of  the 
incoherent  spectrum  referenced  to  the  input.  The  noise  at  sensor  n  is  defined  by 
Gnn(f)  =  [l-Vmn(f)]  Gmm(f)/lHm(f)l2  (20) 

where  Hm(f)  is  the  channel  response  for  sensor  m.  Vmn(f)  is  the  complex  temporal 
coherence  between  sensors  m  and  n 

Vmn(f)  =  [Gmn(f)  Gnm(f)]  /  [Gmm(f)  Gnn(f)]  (21) 

Gmn(f)  is  the  1 -sided  ensemble  average  of  (cross-)  PSD’s  calculated  from  a  window  of 
period  T  by 

Gmn(f)  =  (2/T)  <Xm(f,T)  Xn(f,T)>  (22) 

NY-NEX  channel  hardware  noise  was  measured  between  vertical  sensors  6  and  7 
reinstalled  with  a  0.3  m  separation  at  WO.  The  experiment  recorded  5.12  sec  samples 
every  hour  for  64  hours.  These  samples  were  used  to  construct  the  ensemble  average  of 
periodograms  Gmn(f)  (128  degrees  of  freedom)  graphed  as  the  lower  curves  shown  in  the 
plots  of  Figure  13. 

4.2.2  Ambient  Noise 

Ambient  ground  motion  may  be  defined  as  the  sum  of  independent  random  seismic 
events.  Random  events  may  be  described  as  noise  generated  by  vehicles,  animals,  and  all 
random  electronic  noise  not  described  by  the  channel  response.  Typical  noise  levels 
recorded  during  pre-shot  windows  are  shown  in  Figure  10.  Noise  is  described  both  at  a 
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point  (seismometer)  or  by  that  which  is  passed  by  the  array.  Random  noise  is  classified 
by  its  stationarity,  normality,  and  spectra.  Microseismic  activity  is  found  to  exhibit  the 
attributes  of  a  stationary  Gaussian  random  process. 

Ambient  seismic  noise  has  been  analyzed  for  pre-shot  windows.  In  cases  where 
the  seismic  noise  at  a  point  was  determined  to  be  both  stationary  and  Gaussian,  it  was 
classified  in  terms  of  its  PSD  [16].  Ensemble  averages  of  periodograms  were  calculated 
by  taking  the  sum  of  FFT’s  in  non-overlapping  segments  in  pre-shot  ambient  noise 
windows.  The  data  were  padded  with  zeros  (interval  doubling)  to  reduce  leakage  and 
variance  [17]  generated  by  the  FFT  transformation  (Appendix  A).  The  periodograms 
quantify  the  noise  at  various  frequencies  in  the  0.0  to  50.0  Hz  band  that  may  interfere 
with  the  signal. 

4,2.3  Digital  Noise 

One  important  component  of  digital  noise  which  must  be  avoided  is  aliasing. 
Aliasing  results  from  the  digital  sampling  of  data  below  the  Nyquist  frequency  (twice  the 
highest  frequency)  in  the  analog  waveform.  Sample  rates  below  twice  the  Nyquist 
frequency  will  cause  distortion  in  the  recorded  signal.  The  Protection  Ratio  [18]  identifies 
the  effects  of  the  first-fold  alias  of  the  data  within  the  Nyquist  interval.  It  is  calculated 
by  the  ratio  of  the  amplitudes  of  the  in-band  frequencies  (0.0  to  50.0  Hz)  to  those  of  the 
out-of-band  frequencies  (50.0  to  100.0  Hz)  by 

Hm(f)  =  Hm(f)  /  Hm(2fn-f)  (23) 

where  Hm(f)  is  the  channel  response  for  channel  m  and  fn  is  the  Nyquist  frequency.  The 
high  and  low  gain  Protection  Ratios  of  Figure  14  indicate  that  the  data  recorded  above  30.0 
Hz  is  affected  by  aliasing.  The  Rg  waves  analyzed  in  this  study  are  in  the  1.0  to  5.0  Hz 
band  and  are  minimally  effected  by  aliasing. 
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Another  component  of  digital  noise  results  from  the  digitization  of  analog  signals 
by  the  analog-to-digital  converter  (A/D)  and  is  measured  in  terms  of  quantization  (19]. 
Quantization  noise  has  been  calculated  for  GDAS  and  referenced  to  the  input  by 
Gmm(f)  =  GNN/ 1  Hrn(f)  1 2  (24) 

where  Gmm(f)  is  the  quantization  noise  for  channel  m.  GNN  is  the  quantization  level  Q 
produced  by  the  GDAS  15-bit  A/D  converter  distributed  as  white  noise  over  the  GDAS 
50.0  Hz  Nyquist  interval.  GNN,  shown  in  Equation  (25),  is  measured  in  units  of  volts2/Hz. 
GNN  =  (Q7l2)/50.0  (25) 

4.2.4  Summary 

The  plots  of  Figure  13  show  the  channel  hardware  noise  (lower  curve)  compared 
with  the  ambient  seismic  noise  (upper  curve)  for  channels  1,  3,  7,  and  14.  The  spectra 
are  characteristic  of  the  noise  field  present  throughout  the  NY-NEX  experiment.  Analysis 
of  the  plots  in  the  0.0  to  1.0  Hz  band  indicates  a  sharp  increase  in  the  ambient  noise  from 
1.0  E-10  to  1.0  E-6  (mm/sec )2/Hz.  This  increase  correlates  with  the  increase  in  hardware 
noise  and  may  affect  the  azimuth  and  phase  velocity  estimates.  Figure  13  shows  that  at 
frequencies  greater  than  1.0  Hz,  the  ambient  noise  tends  to  level  off  at  about  1.0  E-ll 
(mm/sec)J/Hz  with  the  exception  of  a  few  noise  peaks  within  the  10-13,  13-16,  18-19,  22- 
26  and  29-30  Hz  bands.  These  noise  peaks  are  of  unknown  origin  and  are,  fortunately, 
outside  of  the  0.6  to  3.3  Hz  Rg  wave  band  of  interest  in  this  study  area  The  consistency 
of  these  spectra  with  those  taken  at  different  ambient  seismic  noise  windows  identifies  the 
noise  as  an  ergodic  random  process.  Overall,  the  spectra  show  that  the  signals  can  be 
detected  accurately  in  the  1.0  to  10.0  Hz  band. 

The  Protection  Ratio  plots  of  Figure  14  indicate,  for  example,  that  samples  recorded 
at  40.0  Hz  will  appear  20  times  greater  than  their  aliased  60.0  Hz  values.  The  Protection 
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Ratio  at  1.0  Hz  is  approximately  1.5  E-4.  Overall,  the  Protection  Ratio  favors  signal 
detection  in  the  0.0  to  40.0  Hz  band. 

A  comparison  of  the  quantization  noise  (Figure  15)  to  the  channel  hardware  noise 
(Figure  13)  indicates  that  the  effect  of  digitization  noise  on  seismic  recordings  is  negligible 
since  it  is  approximately  three  orders  of  magnitude  less  than  the  hardware  noise. 

4.3  ARRAY  NOISE 

Array  noise  is  described  as  noise  passed  by  the  anay.  Sections  4.3.1  and  4.3.2 
describe  this  additive  noise  in  terms  of  spatial  aliasing  and  ambient  noise. 

4.3.1  Spatial  Noise  of  the  Array 

Spatial  aliasing  is  one  component  of  hardware  noise  and  is  controlled  by  sensor 
spacing.  Spatial  aliasing  is  defined  in  terms  of  wavenumber  just  as  sample  aliasing  is 
defined  in  terms  of  frequency.  A  Nyquist  wavenumber  Kn  was  calculated  for  a  sensor 
spacing  D  from  the  following  equations: 

K  =  2rt/£ 

Kn  =  2jx/(D/2) 

K  >  Kn  =  ti/D  (26) 

where  £  is  wavelength.  According  to  Table  3-2  and  Equation  (26),  aliasing  for  the  WO 
NY-NEX  array  would  occur  at  wavenumbers  above  0.1571  rad/m  for  the  minimum  20.0 
m  sensor  spacing  between  channels  7  and  6,  and  above  0.0070  rad/m  for  the  maximum 
448.0  m  sensor  spacing  between  channels  7  and  1.  Spatial  aliasing  was  not  a  noise 
problem  in  this  thesis  since  the  Rg  waves  studied  were  within  a  group  velocity  range  of 
1.9  to  3.4  km/sec,  producing  wavenumbers  within  a  range  of  0.0017  to  0.0009  rad/m. 
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4.3.2  Ambient  Noise  of  the  Array 


Ambient  seismic  noise  attributes  are  given  in  terms  of  frequency  and  wavenumber 
and  are  essential  to  the  optimum  array  processing  of  weak  signals.  Figure  16  shows  the 
ambient  seismic  noise  passed  through  the  array.  The  conventional  frequency-wavenumber 
(F-K)  algorithm  used  to  produce  these  results  is  discussed  in  Section  5.1.1.  The  noise  for 
this  plot  is  block  averaged  over  25  non -overlapping  2.56  sec  segments  (50  degrees  of 
freedom)  in  the  pre-shot  ambient  seismic  noise  window.  Peak  power  (maximum 
P(Kx,Ky,f))  results  for  the  ambient  array  noise  are  listed  in  Table  4-2. 

In  Figure  16,  the  plot  of  P(Kx,Ky,1.76)  for  some  ambient  seismic  noise  indicates 
maximum  power  came  from  a  beam  steered  towards  signals  arriving  from  a  135.0°  azimuth 
with  a  phase  velocity  of  4.69  km/sec  (as  shown  by  the  velocity  circle  on  the  plot).  As 
stated  in  Section  3.3.3,  contour  levels  are  measured  in  db  down  from  the  maximum 
computed  level.  The  plot  of  P(Kx,Ky,1.95)  indicates  a  maximum  power  for  a  beam  steered 
to  signals  arriving  from  a  153.4°  azimuth  with  a  velocity  of  6.58  km/sec.  P(Kx,Ky,2.15), 
and  P(Kx,Ky,2.25)  indicate  a  maximum  power  from  beams  steered  to  receive  signals  from 
the  southeast  and  northeast  at  velocities  around  7.0  km/sec.  Overall,  the  central  main  lobes 
(areas  of  contour  levels  less  than  3  db)  in  Figure  16  indicate  that  coherent  noise  passed  by 
the  array  had  a  variety  of  azimuths  and  high  phase  velocities.  These  velocity  estimates  are 
outside  of  that  expected  for  the  1.9  to  3.4  km/sec  Rg  wave  arrivals.  The  additive  effects 
of  this  ambient  noise  with  the  recorded  signal  will  tend  to  show  up  as  high  phase  velocity 
noise  with  scattered  azimuths. 
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Table  4-2 

AMBIENT 

NOISE  ARRAY  RESULTS 

FREQUENCY 

PEAK  POWER 

PHASE  VELOCITY 

AZIMUTH 

(Hz) 

(mm/sec) 2/Hz 

(km/sec) 

(deq) 

1.76 

0.26 

4.69 

135.00 

1  .  35 

0 . 62 

6.58 

153.44 

2.15 

0 . 15 

7.25 

206.57 

2.25 

0.36 

7.59 

333.43 
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5  RG  WAVE  AZIMUTH  AND  PHASE  VELOCITY  ESTIMATIONS  USING  AN  ARRAY 


5.0  INTRODUCTION 

The  problem  being  studied  here  is  that  of  estimating  Rg  wave  azimuth  and  phase 
velocity  as  a  function  of  frequency  for  the  NY-NEX  shots.  As  mentioned  in  Chapter  1, 
the  Rg  wave  bandwidth  is  limited  to  0.6  to  3.3  Hz.  In  a  simply  layered  half-space, 
azimuth  estimates  should  show  a  minimal  deviation  in  azimuth  from  the  great-circle 
azimuth  for  frequencies  in  the  0.6-3. 3  Hz  range.  However,  it  has  been  found  that  azimuth 
estimates  for  the  NY-NEX  study  were  widely  distributed  about  the  great-circle  azimuth 
within  this  frequency  band.  This  scatter  is  analyzed  to  determine  if  it  is  caused  by  lateral 
heterogeneity  in  the  shallow  crustal  structure  and/or  ambient  noise. 

Classically,  propagating  plane  waves  have  been  characterized  in  terms  of  their 
azimuth  and  phase  velocity  by  a  F-K  power  spectral  density  function  P(Kx.Ky.f)  (see 
Lacoss  [9]  and  Capon  [10]).  This  frequency -domain  synthesis  technique  has  advantages 
over  the  corresponding  time  domain  method  primarily  because  of  considerable  reductions 
in  beam  forming  computation  time  [20],  and  secondarily  the  time  domain  assumption  that 
the  signal  is  identical  across  the  array  [21]  is  not  necessarily  true  for  dispersive  surface 
waves. 

The  azimuth  and  phase  velocity  for  NY-NEX  Rg  waves  were  estimated  using  the 
conventional  F-K  algorithm  discussed  in  Section  5.1.1.  Several  other  conventional  and 
high-resolution  frequency-domain  techniques  were  initially  applied  to  the  Rg  waves  recorded 
by  the  WO  NY-NEX  array  in  an  effort  to  improve  upon  the  azimuth  and  phase  velocity 
estimations.  The  criterion  for  deciding  which  technique  to  implement  was  based  on  finding 
the  method  which  minimized  the  amount  of  scatter  in  azimuth  and  phase  velocity  estimates 
within  the  frequency  band  where  the  Rg  wave  is  present.  The  results  of  the  conventional 
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and  high-resolution  algorithm  applications  developed  by  Capon  et  al.  [  10,20,21]  and  Lacoss 
et  al.  (9,22)  are  summarized  in  the  following  sections. 

1.  Conventional  Methods 

1.  F-K  Analysis  on  Raw  Data 

Out  of  the  several  algorithms  tested.  F-K  analysis  provided  the  smallest  deviation 
in  azimuth  and  phase  velocity  estimates  from  the  expected  values  in  the  0.6-3. 3 
Hz  Rg  wave  band. 

2.  F-K  Analysis  on  Filtered  Raw  Data 

A  0.6-3. 3  Hz  third  order  butterworth  filter  applied  to  the  raw  time-domain  data 
before  F-K  analysis  to  reduce  the  deviation  in  azimuth  and  phase  velocity  estimates 
resulted  in  no  improvements  over  the  results  obtained  from  F-K  analysis  on  raw 
data.  Therefore,  the  additional  amount  of  computation  required  was  not  justified. 

3.  F-K  Analysis  on  Stacked  Data 

Stacking  multiple  shots  occurring  at  the  same  location  (shots  7A,  7B  and  shots 
10B,  10C  for  example)  did  not  reduce  the  deviation  of  azimuth  and  phase  velocity 
estimates  over  those  from  the  single  shet  I'll  analysis  nr.  the  raw  data. 

4.  Wideband  Detection  on  Filtered  Data 

Wideband  detection  is  a  covariance  based  method  that  provides  azimuth  and  phase 
velocity  estimates  for  seismic  signals.  On  raw  data,  this  method  provided  an 
average  azimuth  and  phase  velocity  estimate  for  the  entire  seismic  band,  resulting 
in  deviated  estimates  biased  by  noise  and  other  reflected  and  refracted  signals.  A 
filter  (see  Lacoss  et  al.  (22j)  of  0.6-3. 3  Hz  applied  the  raw  time-domain  data 
reduced  the  deviation  in  the  Rg  wave  azimuth  and  phase  velocity  estimates  from 
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the  expected  values.  These  estimates  were  used  as  a  comparison  to  the 
conventional  F-K  calculated  results. 

11.  High-resolution  Methods 

1.  High-Resolution  (Maximum-Likelihood)  F-K  Analysis  (10] 

For  data  in  the  0.6-3. 3  Hz  Rg  wave  band,  high-resolution  F-K  techniques  resulted 
in  a  greater  deviation  in  azimuth  and  phase  velocity  estimates  compared  to  the 
conventional  F-K  algorithm.  This  scatter  is  probably  caused  by  unstable  PSD 

A  A 

estimates  Pmn(f)  [Pmn(f)  is  discussed  in  Section  5.1.1  ].  This  instability  results  from 

A 

Pmn(f)  being  calculated  from  only  one  ensemble  over  the  data  window.  Capon  [10] 
discusses  the  importance  of  the  block  averaging  for  calculating  a  stable,  inverse 

A 

matrix  Pmn(f)  for  the  high-resolution  method.  He  states,  in  fact,  that  there  should 
be  a  number  of  ensemble  averages  at  least  equal  to  the  number  of  sensors.  In 

A 

analysis  of  the  Rg  wave,  estimates  of  Pmn(f)  averaged  over  L-blocks  were  found 
to  accumulate  noise  due  to  the  existence  of  other  unwanted  signals  in  the  window, 
producing  false  power  maxima  in  the  band. 

2.  High-Resolution  Wideband  Detection  [22] 

High-resolution  wideband  analysis  resulted  in  an  increase  of  azimuth  and  phase 
velocity  estimate  deviations  from  the  values  expected  as  compared  to  the 
conventional  wideband  detection  method.  The  poor  results  of  this  high-resolution 
technique  were  also  attributed  to  the  instability  of  the  inverse  operator  discussed  in 
Section  5.1.1. 
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HI.  Comparisons  of  Methods 


Figure  17  shows  four  typical  examples  of  the  contoured  P(Kx,Ky,f)  values  for  the 
data  processed  from  NY-NEX  shot  5.  Phase  velocity  c(f)  and  azimuth  estimates  are  listed 
at  the  top  of  each  plot.  These  plots  show  the  variations  of  azimuth  and  phase  velocity 
estimates  found  throughout  the  data  processing  for  different  frequencies  in  the  0.6  to  3.3 
Hz  Rg  wave  band.  In  each  plot,  azimuth  and  phase  velocity  estimates  were  determined 
from  the  Kx  and  Ky  grid  locations  where  the  power  spectra  P(Kx,Ky,0  is  at  a  maximum. 
A  dashed  circle  centered  at  the  origin  and  intersecting  this  maximum  indicates  the  phase 
velocity.  Contour  levels  are  shown  in  levels  of  db  down  from  the  normalized  maximum 
(see  Table  3-5  for  description  of  db  levels).  The  azimuth  was  determined  as  the  angle 
between  the  positive  vertical  axis  (north)  and  where  the  phase  velocity  circle  intersects  the 
maximum. 

Figures  18A,  18B,  and  18C  show  results  for  the  wideband  analysis  methods  for 
shot  7 A.  Figure  18B  indicates  a  heightening  of  contours  due  to  the  application  of  a  0.6 
to  3.3  Hz  filter.  The  high-resolution  wideband  results  shown  in  Figure  18C  produced  an 
unacceptable  azimuth  of  around  250°  compared  to  the  expected  301°  great-circle  azimuth. 

Figures  18D,  18E,  and  18F  demonstrate  how  stacking  the  data  effects  F-K  analysis. 
If  noise  was  a  dominant  factor,  azimuth  and  phase  velocity  estimates  would  be  inconsistent 
between  different  recordings  of  data  from  the  same  shot.  Only  minor  differences  are  found 
between  Figures  18D  and  18E  (shown  by  the  arrows  on  the  graphs),  indicating  the 
consistency  in  azimuth  and  phase  velocity  estimates  using  the  F-K  algorithm.  Furthermore. 
Figure  18F  does  not  show  any  improvement  over  Figures  18D  and  18E  from  a  stacking 
of  the  data  from  shots  7A  and  7B  indicating  noise  is  not  a  dominant  factor. 
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5.1  DATA  PROCESSING 


5.1.1  Frequency-Wavenumber  Spectra 

The  conventional  F-K  algorithm  is  used  in  this  thesis  to  estimate  power  as  a 
function  of  frequency  and  wavenumber  P(Kx,Ky,f).  P(Kx,Ky,f)  is  used  to  calculate  the 
azimuth  and  phase  velocity  of  Rg  waves.  It  is  estimated  for  an  NS  sensor  array  by 

NS  A 

P(Kx,Ky,f)  -  (1/NS2)  £  Pmn(u)  exp[jKx(xm-xn)+jKy(ym-yn)]  (27) 
m,n=l 

A 

where  u  is  the  normalized  frequency  u  =  27tk/N,  and  Pmn(u)  is  a  matrix-based  estimate 
of  the  two-dimensional  oss-power  spectral  density  function  between  sensors  m  and  n. 

A 

Pmn(u)  is  calculated  for  L  ensemble  (L-block)  averages  by 
A  L 

Pmn(u)  =  1/L  I  Sml(u)  Snl(u)  (28) 

1=1 

where  Sml(u)  and  Snl(u)  are  the  Fourier  transforms  of  a  N-point  time  series  X  in  the  I'th 
block.  For  sensor  m,  the  discrete  Fourier  transform  (19)  of  Sml(u)  in  the  l’th  block  is 
calculated  by 

N-l 

Sml(u)  =  1/N  I  Xm,i+(1-1)  exp(jiu)  (29) 

i=0 


A 

Pmn(u)  is  assumed  to  be  normalized  to  remove  the  effects  of 
improper  sensor  equalization  by 

A  AAA 

Pmn(u)  =  Pmn(u)/[Pmm(u)  Pnn(u)]l/2  (30) 

Capon  [10J  defines  the  high-resolution  method  for  estimating  the  F-K  spectrum  as 
NS  A 

P’(Kx,Ky,fk)  =  (  E  Qmn(u)  exp[jKx(xm-xn)+jKy(ym-yn)]  )  1/2  (31) 

m,n=l 
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A  A 

where  Qinn(f)  is  the  matrix  inverse  of  the  complex  spectral  matrix  Pmn(f).  The  step  by 
step  procedure  used  in  this  thesis  to  calculate  the  conventional  F-K  spectra  for  the  WO 
NY-NEX  seismic  data  is  summarized  as  follows: 

1.  Determine  the  Group-Velocity  Window  for  Analysis 

The  shot  to  site  distance  R  divided  by  the  maximum  Umax  (3.4  km/sec)  and 
minimum  Umin  (1.9  km/sec)  Rg  wave  group  velocities  determined  the  start  time 
(tO  =  R/Umax)  and  stop  time  (tl  =  RAJ  min)  of  the  group- velocity  window. 

2.  Fourier  Transform  the  Time  Window 

The  FFT  of  each  sensor  was  calculated  for  the  group-velocity  window.  Each  time 
series  was  padded  with  zeroes  (see  Appendix  A)  to  improve  upon  the  estimation 

A 

of  Pmn(u). 

3.  Calculate  the  Cross-Power  Matrix  Pmn(u) 

A 

A  NS  by  NS  (14  by  14)  normalized  cross-power  matrix  of  Pmn(u)  was  generated 
at  frequencies  f  within  the  Rg  wave  band  for  the  calculation  of  P(Kx,Ky,f). 

4.  Calculate  the  F-K  Spectra  P(Kx,Ky,f) 

P(Kx,Ky,f)  was  calculated  over  a  61  by  61  wavenumber  grid  having  a  range  of 
±0.025  rad/m.  The  wavenumber  range  of  ±0.025  rad/m  was  determined  in  order 
to  maintain  accurate  contouring  between  wavenumber  elements  as  shown  in  Figure 
17.  A  61  by  61  grid  partem  for  this  range  results  in  a  phase  velocity  resolution 
of  ±0.4  m/sec. 

5.1.2  Azimuth  and  Phase  Velocity 

The  conventional  F-K  algorithm  estimates  power  P(Kx,Ky,f)  as  a  function  of 
frequency  and  wavenumber,  the  quantities  used  to  calculate  azimuth  and  phase  velocity. 
P(Kx,Ky,f)  was  computed  for  a  Kx  and  Ky  grid  of  wavenumbers  between  ±0.025  rad/m 
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in  the  0.6  to  3.3  Hz  Rg  wave  band.  The  frequency  resolution  within  the  band  is 
determined  by  the  number  of  points  N  in  the  group-velocity  window  (see  step  1  :n  Section 
5.1.1)  and  the  sample  spacing  dt  by 

f  =  K/(Ndt)  (32) 

At  each  frequency  f  within  the  Rg  wave  band,  the  maximum  value  of  P(Kx,Ky,f)  was 
picked  from  the  wavenumber  grid.  The  Kx  and  Ky  location  of  this  maximum  was  used 
to  calculate  the  azimuth  0(f) 

<D(f)  =  arctan(Kx/Ky)  (33) 

and  phase  velocity  c(f) 

c(f)  =  2rrf/lKl  (34) 

where  the  wavenumber  vector  K  was  calculated  from  its  Kx  and  Ky  components  by 

iKl  =  (Kx2  +  Ky2)1"  (35) 


5.1.3  Spatial  Coherence 

The  spatial  coherence  provides  an  analysis  of  coherency  of  the  signals  between 
the  sensors  of  the  WO  NY-NEX  array.  The  spatial  coherence  was  calculated  for 
frequencies  within  the  Rg  wave  band  by 

NS  A 

V(Kx,Ky,f)  =  IP(Kx,Kyd)l  /  I  |Pmn(f)l  (36) 

m,n=l 


5.1.4  Power  Spectral  Density  Sum 

The  power  spectral  density  (PSD)  sum  provides  an  estimate  of  the  average  power 
at  a  frequency  f  for  the  Rg  wave.  The  PSD  sum,  PSD(f),  was  calculated  by  adding  the 
individual  PSD’s  for  each  of  the  NS  sensors  in  the  group-velocity  window  by 
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(37) 


NS  A 

PSD(f)  =  1/NS  I  IPm(f)l 
m=l 


PSD(f)  represents  the  ideal  F-K  power  spectral  density  function.  This  is  shown  in  Equation 

A 

(27)  when  the  exponent  is  zero  (where  Pmn(f)  is  maximum).  Thus,  this  estimate  provides 
information  about  the  amount  of  power  expected  at  a  frequency  f  for  the  Rg  wave  recorded 
by  the  WO  array. 
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6  AZIMUTH  AND  PHASE  VELOCITY  RESULTS 
6.0  INTRODUCTION 

The  NY-NEX  data  interpretation  described  in  this  chapter  was  based  on  the 
combined  analysis  of  (a)  raw  data  records,  (b)  raw  data  records  with  a  0.6-3. 3  Hz  3rd 
order  butterworth  filter  applied,  (c)  ambient  and  hardware  noise  spectra,  (d)  the  PSD  sum 
curve,  and  (e)  plots  of  azimuth  and  phase  velocity  estimates  versus  frequency. 

The  data  were  compared  with  their  average  estimates  and  standard  deviations 
calculated  at  0.5  Hz  intervals  over  a  1.0  to  3.0  Hz  band.  This  procedure  made  it  possible 
to  quantify  the  amount  of  scatter  in  the  azimuth  and  phase  velocity  estimates  at  different 
frequencies  in  the  Rg  wave  band  as  measured  by  the  WO  array.  Theoretically,  the  phase 
velocity  is  not  expected  to  be  less  than  the  group  velocity  at  a  given  period,  and  therefore, 
values  of  phase  velocity  (and  their  associated  azimuths)  are  rejected  if  they  are  less  than 
the  group  velocity.  Spectral  analysis  of  the  NY-NEX  shots  was  limited  to  a  1.0  to  3.0  Hz 
band.  It  was  found  that  for  frequencies  less  than  1.0  Hz  the  data  was  dominated  by- 
additive  noise  (see  Section  4.2.4)  and  for  frequencies  greater  than  3.0  Hz  coherent  Rg 
waves  were  not  present  in  the  spectra. 

NY-NEX  shots  5,  6,  7,  and  8  were  the  only  shots  in  the  study  that  produced 
coherent  Rg  wave  azimuth  and  phase  velocity  estimates  from  the  conventional  F-K  analysis. 
Shot  10  was  processed  to  determine  the  effects  of  stacking  noisy  records  for  improved 
signal-to-noise  ratios  (SNR).  SNR  is  defined  in  this  thesis  as  the  ratio  of  the  signal 
amplitude  specua  lS(f)l  (measured  from  the  spectra  of  the  group-velocity  window)  to  the 
noise  amplitude  spectra  lN(f)l  (measured  from  the  spectra  of  an  ambient  noise  window 
immediately  preceding  the  shot)  as  discussed  in  Section  4.0. 
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6.1  SHOT  ANALYSIS 


6.1.1  Shot  6 

Shot  6,  2000  lbs  in  size,  was  located  at  a  range  of  26.40  km  and  an  azimuth  of 
32.43°  from  the  center  of  the  WO  array.  The  Rg  wave  generated  by  this  shot  appears 
approximately  10.8  sec  after  shot  time  and  can  be  seen  in  both  the  raw  and  filtered  field 
seismograms  (Figures  19A  and  19B).  The  filtered  records  in  Figure  19B  clearly  show  the 
dispersion  of  the  Rg  wave  between  9.8  and  11.8  sec  after  shot  time.  Figure  19C  shows 
a  spectral  comparison  between  the  signal  and  ambient  noise  ground  motions,  and  the 
hardware  noise.  These  spectra  clearly  indicate  a  Rg  wave  in  the  1.0-3.0  Hz  band  with 
amplitudes  ranging  from  0.9  E-7  to  4.0  E-7  (mm/sec)J/Hz  and  a  high  SNR  of  100  to  1. 
The  spectra  have  a  high  spatial  coherence  (0.98  to  1.00)  in  the  1. 0-2.2  Hz  band  as  shown 
in  Figure  19D. 

The  azimuth  and  phase  velocity  estimates  plotted  in  Figure  19D  are  compared  to 
the  azimuth  and  phase  velocity  averages  calculated  at  0.5  Hz  intervals  as  listed  in  Table 
6-1.  Deviation  of  the  azimuth  estimate  was  minimal  (4.01°)  in  the  1. 5-2.0  Hz  band.  The 
high  deviation  (11.34°)  in  the  1.0- 1.5  Hz  band  is  attributed  to  the  high  56.31°  azimuth 
estimate  at  1.17  Hz.  This  scattered  azimuth  corresponds  to  a  decrease  in  PSD(f)  from  2.0 
E-8  to  1.0  E-8  (mm/sec)7Hz  as  shown  in  Figure  19D.  There  is  an  increase  in  azimuth 
scatter  in  the  2.2-2.5  Hz  band,  indicated  by  deviations  of  17.09°.  This  scatter  may  be 
explained  by  a  weak  signal  which  is  shown  as  the  decrease  in  average  power  in  the  PSD(f) 
at  2.5  Hz.  The  phase  velocity  estimates  deviaied  from  their  0.5  Hz  averages  in  the  1.0- 
2.5  Hz  band  by  0.45,  0.30,  and  0.36  km/sec.  These  deviations  are  independent  of  the 
scatter  in  azimuth  estimates,  spatial  coherency,  and  PSD(f)  levels.  Most  of  the  phase 
velocity  estimates  in  the  2.5-3.0  Hz  band  were  rejected  because  they  were  below  the  1.9- 
3.4  km/sec  group-velocity  window.  The  consistently  low  phase  velocities  in  this  band  are 
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associated  with  an  increase  in  signal  amplitude  (indicated  by  PSD(f)  increasing  from  1.0 
E-8  to  6.0  E-8  (km/sec):/Hz)  and  may  be  explained  by  the  presence  of  other  coherent 
waves  in  the  group-velocity  analysis  window. 

An  average  of  the  azimuth  and  phase  velocity  estimates  in  the  1. 0-3.0  Hz  band 
indicates  that  a  beam  was  steered  to  receive  Rg  waves  from  a  49.28°  azimuth  with  a  2.47 
km/sec  phase  velocity.  This  azimuth  is  16.85°  greater  than  the  expected  32.43°  great- 
circle  azimuth.  Figure  1A  indicates  most  of  the  raypath  from  the  shot  point  to  the  array 
as  being  isolated  in  the  2.5  km/sec  Ammonoosuc  Volcanics.  The  French  Pond  granite 
intrusion  is  also  located  between  shot  6  and  the  WO  array  site  and  it  is  thought  to  have 
a  higher  velocity  of  2.7  km/sec  compared  to  the  surrounding  formation.  This  intrusion  may 
have  caused  some  deviation  of  the  Rg  wave  azimuth  from  the  great-circle  path. 
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Table 

6-1 

Shot  6  -  LINEAR  COMPARISONS 

f 

<&(f> 

cD(f )  a 

c  ( f )  c  ( f ) 

o 

(Hz) 

(deg) 

(deg)  (deg) 

(km/sec)  (km/rec 

)  (km/sec) 

1.07 

26.57 

3.62 

1  .17 

56.31 

2.45 

1.27 

33.69 

36.79  11.34 

2.65  2 

.93 

0.45 

1 . 37 

33.69 

2.86 

1  .46 

33.69 

3.06 

1.56 

36.87 

2.36 

1  .  66 

45.00 

2 . 95 

1 .76 

45.00 

42.11  4.01 

2.34  2 

.49 

0.30 

1.86 

38 . 66 

2 .18 

1  .  95 

45.00 

2.60 

2.05 

45.00 

2.73 

2 . 14 

38.66 

2.53 

2.25 

50 .19 

52.83  17.09 

2.17  2  . 

34 

0.36 

2.34 

77 . 47 

1  .  92 

2 .44 

94.76 

1.53 

2.54 

90.00 

1  .  60 

2 . 64 

68.20 

1.85 

2.73 

48.81 

44.05  3.65 

1.94  2. 

01 

0.06 

2.83 

41.19 

2.01 

2.93 

41.19 

2.08 

3.03 

45.00 

2.02 

Average 

azimuth 

<I>(f)  and  phase  velocity  c(f) 

estimates 

are  listed  at  0 

.5  Hz  increments  between  1 

.07 

and  3.03  Hz. 

Anomalous  values  at  2.44-2.64 

Hz  were  omitted 

in  the 

calculations  of 

the  averages 

due  to  their 

low 

phase 

velocity.  G  is 

the  standard  deviation. 
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6.1.2  Shot  7 


Shots  7A  and  7B  were  located  21.33  km  from  the  center  of  the  WO  array  at  an 
azimuth  of  301  /’O0.  "Hie  results  of  the  data  analysis  foi  shots  7A  and  7B  were  nearly 
identical  and  therefore  are  discussed  together.  The  Rg  wave  is  visible  in  the  raw  and 
filtered  time  series  shown  in  Figure  20A  and  20B  at  9.20  sec  after  shot  time,  resulting  in 
a  group  velocity  of  2.32  km/sec.  Examples  of  the  power  spectra  (referenced  to  the  input) 
for  channels  1,  3,  7,  and  14  are  shown  in  the  upper  graphs  of  Figure  20C.  The  Rg  wave 
spectra  appear  in  these  figures  in  the  1. 0-6.5  Hz  band  with  amplitudes  ranging  from  4.0 
E-8  to  3.0  E-5  (mm/sec)2/Hz,  much  higher  than  that  of  shot  6.  The  high  SNR  in  Figure 
20C  (ratio  of  the  upper  curve  to  middle  curve)  is  approximately  100  to  1.  The  noise  in 
the  plots  of  Figure  20C  appears  to  be  uncorrelated  with  the  signal.  Figure  20D  compares 
the  array  analysis  plots.  The  spatial  coherency  is  high  (between  0.97  and  1.00)  in  the  1.1- 
2.2  Hz  band. 

The  azimuth  plot  of  Figure  20D  is  particularly  interesting  because  of  the  changes 
in  azimuth  and  phase  velocity  estimates  with  frequency  in  the  2.0-2.5  Hz  band.  The 
relatively  good  azimuth  and  phase  velocity  estimate  deviations  from  the  average  of  5.84° 
and  0.47  km/sec  in  this  band  correspond  to  a  spatial  coherency  of  0.75  to  0.95.  Shots  7A 
and  7B  contain  similar  power  spectra  characterized  by  erratic  coherence  (8.6  to  9.7)  in  the 
2.25-3.00  Hz  band.  Estimates  in  this  band  were  rejected  on  the  basis  that  their  phase 
velocity  is  below  the  1. 9-3.4  km/sec  group-velocity  window.  The  stack  of  shots  7A  and 
7B  demonstrated  a  slight  improvement  in  azimuth  calculations,  while  the  phase  velocity 
estimates  resulted  in  the  same  amount  of  scatter  (see  Figure  18D,  18E,  and  18F).  This 
comparison  indicates  the  consistency  of  the  F-K  analysis  methods  and  the  minimal  effects 
that  ambient  noise  has  on  these  estimates. 
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Ail  average  of  the  azimuth  and  phase  velocity  estimates  in  the  1. 0-2.5  Hz  band 
indicates  that  a  beam  was  steered  to  receive  the  Rg  waves  from  a  276.74°  azimuth  for  a 
2.62  km/sec  phase  velocity.  This  azimuth  is  24.46°  less  than  the  expected  great-circle 
azimuth.  This  change  in  azimuth  could  have  been  caused  by  a  rock  velocity  gradient 
decreasing  at  an  angle  to  the  direction  of  propagation  of  the  Rg  wavefront.  A  simple 
model  to  demonstrate  this  effect  is  to  consider  the  propagation  of  Rg  waves  across  a 
boundary  between  two  different  seismic  velocities.  For  a  high  velocity  (VI)  to  low 
velocity  (V2)  boundary,  the  angle  of  incidence  il  is  greater  than  the  angle  of  refraction  i2 
(see  Figure  21).  The  24.46°  raypath  deflection  is  consistent  with  Snell’s  law  in  Equation 
(38)  for  a  high-to-low  velocity  boundary  somewhere  along  the  raypath 

V1/V2  =  sin(il)/sin(i2)  (38) 

The  geologic  structure  between  shot  7  (West  Groton,Vt)  and  the  array  is  best  shown  by 
the  cross-section  in  Figure  IB.  This  cross-sections  shows  several  Ordovician  sequences 
(with  velocities  ranging  from  2.5  to  3.0  km/sec)  contacting  with  the  lower  velocity  (2.5 
km/sec)  Ammonoosuc  Volcanics  to  the  southeast.  This  contrasting  geology  could  explain 
the  deviation  from  the  great-circle  raypath  seen  in  the  observed  data. 
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Table 

6-2 

Shot 

7  -  LINEAR  COMPARISONS 

f 

6(f) 

6(f) 

a 

c  ( f ) 

c  (f ) 

a 

(Hz) 

(deq) 

(den) 

(deq) 

(km/sec) 

(km/sec) 

(km/ sec ) 

1  .  07 

90 . 00 

4.05 

1 .17 

270.00 

1.77 

1.27 

303.69 

292.50 

15.83 

2 . 65 

2.41 

0  .  34 

1.37 

296.56 

4 . 61 

1.46 

281.31 

2.17 

1.56 

291 . 80 

2.19 

1  .  66 

284 .04 

3.04 

1 .76 

281.31 

283.82 

5.84 

2.60 

2.45 

0.47 

1 .86 

278 .13 

1  .  98 

1  .  95 

270 . 00 

1  .84 

2.05 

251.56 

2.44 

2 .14 

248.20 

3.01 

2.25 

258.69 

264.77 

16.96 

3.32 

2.84 

0.36 

2 . 34 

279.46 

2 . 91 

2 .44 

285.95 

2.53 

2.54 

270.00 

1 . 60 

2.64 

255.96 

1.61 

2 .73 

251.57 

rejected 

1.63 

rejected 

2.83 

262.41 

1 .41 

2 . 93 

265.91 

1.57 

3.03 

266.42 

1  .  42 

Azimuth  <t>  ( f )  and  phase  velocity  c(f)  estimates  are  listed 
at  increments  between  1.07  and  3.03  Hz.  Anomalous  values 
at  1.07,  1.17,  1.37,  1.95,  and  2.54-3.03  Hz  were  rejected 
due  to  their  phase  velocities  being  outside  of  the  group 
velocity  window. 
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6-1-3  Shot  5 


Shot  5,  2000  lbs  in  size,  was  located  57.45  km  from  the  center  of  the  WO  array 
at  an  azimuth  of  60.01°.  Recognition  of  a  coherent  Rg  wave  is  limited  to  channels  5.  6, 
7.  and  10  in  the  raw  and  filtered  time  series  shown  in  Figures  22A  and  22B,  indicating 
spatial  coherency  for  channels  within  a  20  m  separation.  Comparisons  of  the  spectra  for 
channels  1,  2,  4,  and  5  in  Figure  22C  indicate  a  SNR  [ratio  of  the  upper  (signal)  curve  to 
the  middle  (ambient  noise)  curve]  of  20  to  1.  In  the  PSD  sum  curve  of  Figure  22D,  the 
Rg  wave  spectra  is  apparent  in  the  1. 0-2.3  Hz  band.  The  spectra  are  characterized  by  a 
high  spatial  coherence  in  the  range  of  0.92  to  1.00  in  the  1. 0-2.0  Hz  band  and  a  rapid  loss 
of  spatial  coherency  from  0.90  to  0.55  in  the  2.0-3.0  Hz  band. 

The  azimuth  and  phase  velocity  estimates  in  the  1. 0-3.0  Hz  band  plotted  in  Figure 
22D  have  been  compared  with  their  average  estimates  and  standard  deviations  at  0.5  Hz 
intervals  (see  Table  6-3)  to  describe  the  amount  of  scatter.  The  azimuth  plot  of  Figure 
22D  and  the  azimuth  standard  deviations  of  28.55°  indicate  that  the  scatter  is  high  in  the 
2  0-2.5  Hz  band.  Azimuth  estimates  show  a  minimal  deviation  in  the  1 .0-1.5  Hz  band  with 
a  deviation  of  2.76°.  The  scatter  in  phase  velocity  for  the  1. 0-3.0  Hz  band  is  below  0.38 
km/sec.  However,  numerous  phase  velocities  were  rejected  because  they  were  below  1.9 
km/sec  in  the  2.0-3.0  Hz  band.  The  phase  velocity  plot  of  Figure  22D  and  the  average 
estimates  at  0.5  Hz  in  Table  6-3  resemble  a  dispersive  phase  velocity  curve. 

An  overall  average  of  the  azimuth  and  phase  velocity  estimates  in  the  1. 0-3.0  Hz 
band  (avoiding  the  anomalous  values)  indicates  that  a  beam  was  steered  to  receive  Rg 
waves  from  an  azimuth  of  64.63°  with  a  phase  velocity  of  2.60  km/sec.  This  azimuth 
deviation  of  4.62°  indicates  essentially  no  changes  in  the  lateral  variations  in  the  rock 
velocities  along  the  great-circle  azimuth.  Geologically,  shot  5  parallels  the  north-northeast 
trending  structure  in  the  region  (see  Figure  1A). 
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Table 

6-3 

Shot 

5  -  LINEAR  COMPARISONS 

o 

(km/ sec) 

f 

(Hz) 

<D(f) 

(deg) 

<l>(f ) 

(deg) 

CT 

(deg) 

c  (f ) 

(km/sec) 

c  (f ) 

(km/sec) 

1.07 

71 . 57 

2 .56 

1  . 17 

71.57 

2.79 

1  .27 

71  .57 

71.77 

2.76 

3.03 

2.59 

0 .42 

1  .  36 

68.20 

1  .  91 

1.46 

75 . 96 

2 . 68 

1.56 

90.00 

1.47 

1  .  66 

56 . 31 

3.47 

1.76 

63.44 

62.65 

5.98 

2 . 96 

3.05 

0  .  38 

1.86 

108 .44 

2.21 

1  .  95 

68.20 

2.73 

2.05 

78 . 69 

1  .  52 

2 . 15 

33.69 

2.25 

2.25 

74.06 

53 .88 

28.55 

2 . 33 

2.29 

0.06 

2 . 34 

235.30 

1.12 

2 .44 

143.13 

1  .  84 

2.54 

250.91 

0.70 

2 . 64 

209.05 

1  .  93 

2.73 

56 . 31 

56 .31 

0.00 

1  .  91 

1.91 

0.00 

2.83 

81.25 

1  .  62 

2.93 

45.00 

1  .  30 

3.03 

225.00 

0.54 

Average  azimuth  <1>  (f )  and  phase  velocity  c(f)  estimates 
are  listed  at  0.5  Hz  increments  between  1.07  and  3.03  Hz. 
Anomalous  values  at  1.56,  2.05,  2.34-2.54  and  2.83-3.03 
Hz  were  omitted  in  the  calculations  of  the  averages 
because  of  their  low  phase  velocity.  Anomalous  values  at 
1.86  and  2.64  Hz  were  omitted  in  the  calculations  of  the 
averages  because  of  their  low  azimuth. 
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6.1.4  Shot  8 


Shot  8  was  2000  lbs  in  size  with  an  azimuth  of  280.15°  from  the  WO  array  center. 
Shot  8  had  nearly  the  same  azimuth  as  shot  7  (301.20°)  but  had  a  range  of  46.19  km 
instead  of  21.33  km.  Shot  8  shows  evidence  of  the  Rg  wave  in  the  filtered  time  plot  of 
Figure  23B.  The  Rg  wave  appears  at  18.6  sec  from  shot  time  indicating  a  group  velocity 
of  2.6  km/sec.  The  PSD  plots  in  Figure  23C  (upper  curve)  show  an  appearance  of  the  Rg 
wave  in  the  1. 0-3.0  Hz  band  with  a  SNR  of  32  to  1.  The  spatial  coherency  of  these  shots 
is  high  (0.98  to  1.00)  in  the  1. 0-2.0  Hz  band,  decreasing  to  around  0.80  in  the  2. 0-3.0  Hz. 
band.  The  average  power  in  Figure  23D  in  the  1. 0-2.0  Hz  band  is  between  0.8  E-8  and 
1.0  E-8  (mm/sec)2/Hz. 

Azimuth  estimates  in  the  1. 0-2.0  Hz  band  show  a  deviation  of  2.71°  to  9.20°  from 
the  average.  The  phase  velocity  estimates  in  the  1. 0-2.0  Hz  band  deviate  from  their 
averages  by  0.53  to  0.47  km/sec.  This  large  scatter  in  phase  velocity  estimates  is  clear  in 
the  phase  velocity  plot  of  Figure  23D.  Most  of  the  estimates  in  the  2.0-3.0  Hz  band  were 
unacceptable  because  of  their  low  phase  velocities. 

The  average  azimuth  and  phase  velocity  estimates  for  the  1. 0-3.0  Hz  band  indicate 
that  a  beam  was  steered  at  an  azimuth  of  266.62°  to  receive  Rg  waves  having  a  phase 
velocity  of  2.60  km/sec.  This  azimuth  is  34.58°  less  than  the  expected  280.15°  great- 
circle  azimuth.  Such  a  variation  in  the  observed  azimuth  from  the  expected  great-circle 
azimuth  could  have  been  caused  by  lateral  changes  in  rock  velocities  along  the  raypath 
similar  to  that  suggested  to  explain  the  observations  from  shot  point  7. 
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Anomalous  values  at  2.14-2.64,  2.83,  and  2.93  Hz  were 
omitted  in  the  calculations  of  the  averages  because  of 
their  low  phase  velocity.  The  anomalous  value  at  1.27  Hz 
was  omitted  in  the  calculation  of  the  1.0-1. 5  Hz  average 
because  of  its  low  azimuth. 
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6.1.5  Shot  10 


Shot  10  was  located  at  an  azimuth  of  269.84°  and  a  range  of  1 10.47  km  from  the 
WO  array  site.  This  shot  was  near  the  farthest  limit  from  which  Rg  waves  with  acceptable 
amplitudes  for  analysis  were  detected  by  the  array  during  the  N'Y-NEX  experiment.  The 
processing  of  this  shot  was  to  determine  if  the  stacking  of  two  noisy  data  records  from  two 
common  shots  (10B  and  10C)  would  increase  the  SNR,  causing  a  decrease  in  the  scatter 
of  azimuth  and  phase  velocity  estimates  in  the  1. 0-3.0  Hz  band.  The  Rg  wave  is  not 
visually  present  in  the  raw  and  filtered  field  data  records  shown  in  the  group-velocity 
windows  of  Figures  24A  and  24B.  The  plots  in  Figure  24C  also  appear  to  show  no 
presence  of  the  Rg  wave  in  the  1. 0-3.0  Hz  band.  The  low  SNR  in  this  band  is  estimated 
at  2  to  1.  Figure  24D  shows  the  F-K  results  due  to  the  time  series  stacking.  Azimuth 
estimates  in  Figure  24D  have  little  indication  of  a  trend. 

An  average  of  the  azimuth  and  phase  velocity  estimates  over  the  1. 0-3.0  Hz  band 
indicates  that  a  beam  was  steered  to  receive  Rg  waves  from  an  azimuth  of  239.81°  with 
a  phase  velocity  of  2.88  km/sec.  This  azimuth  is  at  a  decrease  of  30.03°  from  the 
expected  azimuth,  indicating  a  shot-to-site  lateral  decrease  in  rock  velocity  along  the  great- 
circle  azimuth.  The  consistency  in  F-K  estimates  between  the  individual  shots  and  their 
stack  indicates  that  ambient  noise  did  not  significantly  effect  the  estimates  of  azimuth  or 
phase  velocity. 
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6.1.6  Summary 


Figures  19  thru  24  represent  the  time  and  spatial  processing  results  for  the  NY- 
NEX  shot  data  recorded  at  the  WO  array.  It  was  found  that  the  stacking  of  two  shots 
from  common  shot  points  did  not  improve  the  F-K  azimuth  and  phase  velocity  estimates, 
suggesting  that  the  scatter  in  these  estimates  was  not  due  to  additive  noise.  Table  6-5 
indicates  the  effect  that  azimuth  measurement  error  *4>rms  and  phase  velocity  measurement 
error  *crms/c  has  on  scatter  (as  discussed  in  Section  4.0).  The  high  RMS  levels  in  phase 
velocity  measurement  error  listed  for  shot  10  may  explain  the  large  amount  of  scatter  in 
the  conventional  F-K  estimates  in  the  1.0  to  3.0  Hz  band.  Table  6-5  indicates  a  minimal 
deviation  in  azimuth  measurement  error  for  all  the  shots  processed. 


Table  6- 

5 

ESTIMATION  ERROR 

SHOT 

SNR 

MJrms 

RANGE  OF 
£/*x 

acrms/c 

5 

20 

0.04 

4.2  -  170.0 

0.03  - 

1 . 35 

6 

100 

0.01 

4.2  -  170.0 

0.01  - 

0.27 

7 

100 

0.01 

4.2  -  170.0 

0.01  - 

0.27 

8 

32 

0.02 

4.2  -  170.0 

0.02  - 

0.85 

10 

2 

0.35 

4.2  -  170.0 

0.33  - 

13.53 

Rg  wavelengths 

+  sensor 

spacings  (3400  m/sec 

+  20  m) 

and  1900  m/sec 
minimum  ranges 

+  448  m) 
of  £/Ax. 

provides  the 

maximum 

and 

For  the  NY-NEX  shot  data  analyzed,  the  1.0  to  2.0  Hz  band  consistently  presented 
a  high  spatial  coherence  (ranging  between  0.95  and  1.00).  The  spatial  coherence  in  the  2.0 
to  3.0  Hz  band  was  found  to  be  below  0.90.  Dispersive  phase  velocity  trends  were  noticed 
in  the  0.5  Hz  incremental  averages  for  shots  5,  6,  7,  and  8  in  Tables  6-1  to  6-4. 

As  a  comparison  to  the  average  F-K  estimates,  the  conventional  algorithm  for 
wideband  analysis  (see  Section  5.0,  1-4)  with  a  0.6  to  3.3  Hz  filter  was  applied  in  the 
group- velocity  window  over  the  0.6  to  3.3  Hz  frequency  band.  Azimuth  and  phase  velocity 
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estimates  from  the  F-K  analysis  and  the  wideband  analysis  techniques  are  compared  in 
Table  6-6,  and  the  wideband  F-K  estimates  are  shown  in  Figure  25.  Wideband  azimuth 
estimates  in  Table  6-6  are  within  10°  of  the  F-K  azimuth  estimates.  The  wideband  phase 
velocity  estimates  in  Table  6-6  are  unacceptable  for  shots  5  and  7B  because  of  their  low 
phase  velocity.  The  wideband  phase  velocity  estimates  for  shots  6  and  7A  are  lower  than 
the  2.5  to  3.0  km/sec  Rg  wave  velocities  expected  in  the  region.  Overall  the  wideband 


estimates  indicate  that  the  phase  velocity  estimates  are  poorly  determined  in  the  Rg  wave 
analysis. 


Table  6-6  AZIMUTHS  AND  PHASE  VELOCITIES 
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A  map  showing  true  azimuths  and  F-K  azimuths  for  NY-NEX  shots  5,  6,  7,  X.  and 
10  is  shown  in  Figure  25.  The  deflection  of  the  Rg  raypaths  shows  a  consistent  pattern 
which  may  reflect  velocity  variations  in  the  local  geology.  A  model  to  explain  the  raypath 
deviations  is  one  with  a  gradient  in  the  rock  velocity,  with  the  velocity  decreasing  to  the 
southeast.  The  Rg  wave  recorded  from  shot  5  was  less  likely  to  be  scattered  by  geologic 
features  such  as  faults  since  its  raypath  is  parallel  to  the  north -northeast  trending  trust  faults 
and  parallel  to  the  strike  of  the  lateral  velocity  change.  Those  raypaths  at  greater  angles 


to  this  northeast  direction  of  shot  5  show  increasing  raypath  deviations  from  the  great- 
circle  path. 

A  decrease  in  rock  velocity  to  the  southeast  between  the  sites  of  those  NY-NEX 
shots  processed  in  this  thesis  and  the  WO  array  site  can  be  attributed  to  a  transition  from 
the  higher  velocity  metamorphic  formations  of  Ow,  Og,  and  Oal  to  the  lower  velocity  Oam 
formation  from  northwest  to  southeast.  F-K  analysis  results  for  shot  6  found  an  increase 
in  the  Rg  raypath  azimuth  of  16.85°.  F-K  analysis  for  shot  5  found  a  small  increase  in 
azimuth  of  4.62°.  From  the  northwest,  shot  7  was  found  through  F-K  analysis  to  have 
decreased  by  24.46°  in  azimuth  from  the  true  azimuth.  The  average  F-K  azimuth  estimate 
for  shot  8  was  266.62°,  34.58°  less  than  the  great-circle  azimuth.  An  average  of  the 
azimuth  and  phase  velocity  estimates  for  shot  10  indicated  a  decrease  in  azimuth  of  30.03° 
from  the  great-circle  azimuth.  These  raypath  deflections  support  a  model  where  shots  5 
and  6  are  not  deflected  much  by  lateral  rock  velocity  variations  since  their  raypaths  were 
parallel  to  the  grain  of  the  structure.  Rg  waves  from  shots  7,  8,  and  10  had  raypaths  at 
different  angles  to  the  geologic  structure,  and  in  all  cases  slower  velocities  towards  the 
southeast  would  have  deflected  the  rayp?'>  in  a  counterclockwise  direction.  This  model 
is  illustrated  in  Figure  25. 

6.2  COMPARISON  OF  GROUP  AND  PHASE  VELOCITY 

Kafka  (personal  communication,  1989)  provided  several  group  velocity  dispersion 
curves  calculated  from  the  data  recorded  at  sensors  1-7  for  shot  22.  These  calculations 
were  based  on  an  algorithm  developed  by  Dziewonski  et  al.  [23]  for  the  analysis  of 
transient  seismic  signals.  The  data  were  processed  at  WO  using  a  narrow-band  pass  filter 
over  a  group-velocity  window  of  2.2  to  3.2  km/sec.  The  results  for  each  sensor  are  shown 
as  the  lower  curves  in  Figure  27.  Figure  27  is  a  typical  example  of  a  Rg  wave  group 
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velocity  dispersion  produced  by  quarry  blasts  in  the  New  England  region  [2J.  The 
difference  between  group  velocity  measurements  for  each  sensor  is  within  0.11  km/sec, 
indicating  the  reliability  of  the  estimates.  Acceptable  phase  velocity  estimates  in  Tables 
6-1  to  6-4  were  averaged  at  0.5  Hz  increments  to  provide  one  general  phase  velocity 
dispersion  curve  c(av)  for  the  Rg  waves  recorded  by  the  NY-NEX  shots.  This  phase 
velocity  dispersion  curve  is  shown  in  Figure  26.  From  c(av),  the  group  velocities  listed 
in  Table  6-7  were  calculated  using  Equation  (2),  Equation  (3)  and  the  derivatives  of  K=w/c, 
resulting  in  the  equation 

dK/dw  =  1/c  -  w/c2  dc/dw  =  1  AT  (39) 

Rewriting  Equation  (39)  as  a  function  of  frequency  f  results  in 

dK/dw  =  l/c(f)  -  f/c(f)2  [c(f)-c(f-df)]/df  =  1/U(f)  (40) 

The  group  velocity  calculations  U’(f)  agree  within  0.3  km/sec  of  the  group  velocity 
dispersion  curve  in  the  1.00  to  2.25  Hz  frequency  band.  U’(2.75)  is  lower  than  the 
expected  U(2.75)  as  listed  in  Table  6-7  and  shown  in  Figure  27.  This  low  value  may  be 
caused  by  the  poor  Rg  signal  generally  found  for  frequencies  above  2.00  Hz.  The  results 
verify  that  Rg  wave  phase  velocity  dispersion  is  present  and  that  the  expected  results  agree 
with  other  New  England  Rg  wave  dispersive  group  velocity  models  produced  by  Kafka  [2]. 


Table  6-7 

GROUP  VELOCITY  MEASUREMENTS 

f 

C  (f ) 

df 

c (f-df) 

U'  (f) 

U(f) 

(Hz) 

(km/sec) 

(Hz) 

(km/sec) 

(km/sec) 

(km/sec) 

1.25 

2.74 

0.25 

2 . 80 

2.47 

2.59 

1.75 

2.57 

0.50 

2.74 

2.09 

2.41 

2.25 

2.56 

0.50 

2.57 

2.52 

2.36 

2.75 

2 .01 

0.50 

2.56 

0 . 82 

2.30 

df  is  the 

frequency  step  between 

c  (f )  and 

c (f-df) , 

U'  (f)  is 

the  group 

velocity  calculated  from  values 

in  this  table  and 

Equation  (40), 

and  U(f) 

is  the 

group  velocity  shown  in 

Figure  27  . 
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7  CONCLUSION 


A  conventional  F-K  algorithm  has  been  used  to  estimate  azimuths  and  phase 
velocities  for  Rg  waves  produced  from  the  NY-NEX  experiment.  Data  analysis  was  based 
on  producing  estimates  of  P(Kx,Ky,f)  over  a  1.0  to  3.0  Hz  band,  which  was  shown  by 
Kafka  [2]  to  contain  Rg  waves  in  New  England.  The  noise  at  the  WO  NY-NEX  array  was 
identified  in  an  effort  to  determine  if  the  scattered  azimuth  and  phase  velocity  estimates 
were  the  result  of  (a)  scattering  due  to  lateral  heterogeneity  in  geology  or  (b)  measurement 
error  due  to  the  ambient  seismic  noise. 

The  scatter  in  azimuth  and  phase  velocity  estimates  appears  to  be  primarily  due  to 
a  lateral  variation  in  the  shallow  crust  and  not  noise.  This  lateral  variation  appears  to  be 
related  to  variations  in  the  seismic  velocity  structure  of  the  rocks  underlying  the  various 
geologic  formations  outcropping  in  the  study  area.  These  conclusions  are  based  on  (a)  the 
RMS  error  in  phase  velocity  measurements  between  two  sensors  was  found  to  be  negligible 
for  the  near  shots  (5,  6,  7,  and  8),  but  significant  for  the  far  shot  (10),  (b)  SNR 
improvements  resulting  from  the  stacking  of  common  data  from  two  shots  at  the  same 
shot  point  was  shown  to  not  reduce  the  scatter  in  the  azimuth  and  phase  velocity  estimates 
at  different  frequencies,  (c)  ambient  noise  was  identified  through  conventional  F-K  analysis 
to  be  multidirectional  with  phase  velocities  in  the  4.7  to  7.6  km/sec  range,  having  no 
additive  effects  on  the  Rg  wave,  and  (d)  Rg  raypath  azimuth  estimates  agree  with  a  model 
having  a  laterally  decreasing  rock  velocities  to  the  southeast. 

Azimuth  and  phase  velocity  estimates  in  the  Rg  wave  1.0  to  3.0  Hz  band  were 
found  to  be  affected  by  the  regional  geology.  Shot  5,  with  a  raypath  parallel  to  the  north- 
northeast  trending  structure,  was  found  to  be  deflected  by  only  5°  over  the  great-circle 
azimuth.  The  deflections  found  in  shots  6,  7,  and  8  are  characteristic  of  lateral  velocity 
variations  along  the  great-circle  azimuths.  F-K  analysis  resulted  in  average  phase  velocity 
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estimates  around  2.6  km/sec,  within  the  2.5  to  3.0  km/sec  group  velocity  range  expected 
for  the  metamorphic  to  igneous  rocks  of  the  region. 
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APPENDIX  A:  THE  POWER  SPECTRAL  DENSITY  FUNCTION 


The  Power  Spectral  Density  (G(f)  or  PSD)  is  a  2-sided  even  function  with  one  half 
the  total  power  (variance  per  cycle  per  second)  in  the  positive  frequency  range.  In  this 
paper,  the  PSD  considered  over  positive  frequency  is  the  one-sided  PSD  function  2S(f). 
Thus,  the  PSD  is  expressed  in  terms  of  a  density  2S(f)  for  positive  frequencies. 

To  expand  on  this  matter  theoretically,  S(f)  is  defined  in  terms  of  the 
autocorrelation  function  R(x) 


SCO  =  /  R(t)  exp(-jwx)  di  (A-l) 


and  relating  to  the  time  series  x(t) 

T 

S(f)  =  lim  1/(2T)  I  |  x(t)  exp(-jwt)  dt  1 1  (A-2) 

T-*»  -T 

S(f)df  describes  the  contribution  to  the  variance  (amount  of  power)  within  the  band  f  and 
f+df.  If  x(t)  is  real,  as  in  the  case  of  the  time  history  records,  S(f)  may  be  represented  by 
a  2-sided  even  function 

T 

S(f)  =  lim  (1/2T)  I  j  x(t)  cos(wt)  dt  l:  (A-3) 

T-**>  -T 

The  one-sided  PSD  is  shown  as 


T 


G(f) 

=  lim  (1/T) 

1  {  x(t)  cos(wt)  dt  iz 

(A-4) 

T— 

0 

G(0 

=  2  S(f) 

f  >  0 

=  0 

f  <  0 

(A-5) 

The  PSD  is  defined  over  the  interval  [-T,+T],  and  the  Fourier  transform  of  x(t)  is 
defined  over  the  interval  [0,+T].  To  obtain  the  true  PSD  representation  of  the 
autocorrelation,  the  interval  of  x(t)  is  doubled.  To  accomplish  this  task.  x(t)  is  increased 
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APPENDIX  A:  THE  POWER  SPECTRAL  DENSITY  FUNCTION 


from  [0,+T]  to  [0,2T]  with  the  installation  of  T  zero  pads.  Finally,  Equation  (A-4)  is 
written  in  terms  of  the  autocorrelation 

OO 

G(f)  =  2  /  R(t)  cos(wt)  dx  (A-6) 

0 

Relating  the  above  one-sided  power  spectra  to  the  two-sided  power  spectra,  the  total  power 
in  the  positive  frequency  range  contains  the  power  in  the  negative  frequency  range  plus  the 
power  in  the  positive  frequency  range. 
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APPENDIX  B:  SYMBOLS 


Variable  Description 

p  Mean  (Estimated  Mean  of  the  Sample) 

u  Normalized  Frequency  u  =  2nrc/(N-l) 

a  Standard  Deviation 

cr  Variance 

<&(f)  Azimuth  as  a  Function  of  Frequency 


71 

<> 

l/(Ndt) 

a 

A(z) 

B(Kx,Ky,t) 

c(0 

D 

df 

dt 

dt 

dw 

E 

E 

exp 

f 

fn 

F 


Pi  (7T  =  3.14159  26536) 

Expected 

Fundamental  Frequency 
Acceleration 

Amplitude  as  a  Function  of  Depth 

Spatial  Window  Function  or  Beamforming  Array  Response 
Phase  Velocity 
Sensor  Spacing 

Frequency  Spacing  between  Samples 
Time  Spacing  between  Samples 
Variable  Spacing 

Wavenumber  Spacing  between  Samples 
Base  10  Exponent  (E-2  is  equivalent  to  102) 

Frequency  Band 

Exponential  Function  (exp  =  2.71828  18285) 

Frequency,  Discrete  f  =  K/(Ndt) 

Nyquist  Frequency 
Force 
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APPENDIX  B:  SYMBOLS 


g 

Gmn(f) 

Gnn(f) 

GNN(f) 

H(0 


Gravity 

1 -Sided  (Cross-)  Power  Spectral  Density 
Incoherent  Noise 
Hardware  Noise 

Channel  Response  [(OUTPUT)/(INPUT)  =  (volts)/(mm/sec)] 
Current 


j 

K 

km 

Kn 

L 

m 

M 

n 

N 

N(f) 

NS 

P(f) 

P(Kx,Ky,f) 

P(x) 

P(x) 

Pmn(u) 

PSD 


ft 

Wavenumber  (rad/m) 

Kilometers 

Nyquist  Wavenumber  (;t/D) 

Number  of  Time  Segment  Elocks  (Windows) 

Meters,  or  Sensor  Channel  Number 
Mass 

Sensor  Channel  Number 

Number  of  Samples  in  Time  Window 

Noise  Spectra 

Number  of  Sensors 

1 -Sided  Power  Spectral  Density  (PSD),  G(f) 

Frequency-wavenumber  Power  Spectral  Density 

Probability  Density  Function 

Probability  Distribution  Function,  Prob(x  <  x’) 

Spectral  Matrix  Calculated  from  the  Discrete  Fourier  Transform  of 
the  Cross-correlation  Function 
Power  Spectra]  Density,  P(f) 
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Q  Quality  Factor 

A  A 

Qmn(u)  Inverse  of  the  Spectral  Matrix  Pmn(u) 

R(t)  Autocorrelation  Function 

S(f)  2-Sided  Power  Spectral  Density 

T  Period  (T  =  I/O 

U  Group  Velocity 

u(x,h,z,t)  Panicle  Displacement  in  Canesian  Coordinates 

V(Kx,Ky,0  Spatial  Coherence  Function 

Vmn(0  Temporal  Coherence  Function 

Vp  P-Wave  Velocity 

Vr  Rayleigh  Wave  Velocity 

Vs  S-Wave  Velocity 

w  Angular  Frequency  (rad/sec) 

w(x,y),Wmn  Array  Weighting  Function 
x  Random  Variable  or  Cartesian  Coordinate 

xm  x-coordinate  Position  of  Sensor  m 

xn  x-coordinate  Position  of  Sensor  n 

x’  Random  Variate 

x(t)  Time  Series 

y  Cartesian  Coordinate 

ym  y-coordinate  Position  of  Sensor  m 

yn  y-coordinate  Position  of  Sensor  n 

l  Cartesian  Coordinate 
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The  diagram  shows  the  Geophysical  Data  Acquisition  System  (GDAS)  hardware  flow 
chart.  The  "signal  conditioner"  module  contains  the  preamplifier  and  filter  electronics. 


Figure  7 
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UlPfcgB  SHOT  18  UU.9-3.4i  *HS£C  3gO  BUTFEWWOE  TH  PftSS  I  3  H7 _  Ulpt>*»  SHOT  iB  nil  .9-3.41  Kf1  St  I  3RQ  BUT  ItfiWUjR  1 H  PmSS 
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OOLTS 


C^uitw J  FREQUENCY  DONAiH  PLOTTING  GfUlEU:  FREQUENCY"  UJfiaiN  PlOTTINC. 


FIGURE  23C 
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»PfcQUfcl*Y  IK.’ 


4*>  .  i(*  fc  H  K*’  inUFH  t’Htl  .  1  *■» 


FIGURE  24B 
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(?Uf£Mi  FfitQUtMCY  DOMAIN  PLOTTING  ffcfOUfNTV  fuTNA  1 M  Pi  CiT  TING 


FIGURE  24C 
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FIGURE  24D 
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